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Numerical  Modeling  of 
Micromechanical  Devices 
Using  the  Direct  Simulation 
Monte  Carlo  Method 

A  direct  simulation  Monte  Carlo  ( DSMC)  investigation  of  flows  related  to  microelec¬ 
tromechanical  systems  (MEMS)  is  detailed.  This  effort  is  intended  to  provide  tools 
to  facilitate  the  design  and  optimization  of  micro-devices  as  well  as  to  probe  the 
effects  of  rarefaction,  especially  in  regimes  not  amenable  to  other  means  of  analysis. 
The  code  written  for  this  purpose  employs  an  unstructured  grid,  a  trajectory-tracing 
particle  movement  scheme,  and  an  infinite  channel’*  boundary  formulation.  Its 
results  for  slip-flow  and  transition  regime  micro-channels  and  a  micro-nozzle  are 
presented  to  demonstrate  its  capabilities. 


Introduction 

Microelectromechanical  systems  (MEMS)  are  the  subject  of 
increasingly  active  research  in  a  widening  field  of  disciplines. 
These  are  veiy  small  (micron-scale)  sensors  and  actuators  man¬ 
ufactured  with  photolithographic  techniques  similar  to  those 
used  for  integrated  circuit  (IC)  chips.  MEMS  applications  range 
from  consumer  products  (air-bag  triggers,  micro-mirror  dis¬ 
plays),  to  industrial  and  medical  tools  (micro- valves,  micro¬ 
motors),  to  instrumentation  (micro  pressure  sensors,  micro 
shear-stress  sensors)  (Scott,  1993). 

The  small  size  of  MEMS  poses  unique  challenges  in  the 
design  phase,  however.  While  the  mechanical  properties  of  mi- 
cromachined  materials  are  reasonably  well-studied,  fluid  effects 
at  micron  scales  are  not.  These  effects,  such  as  film  damping 
of  resonant  structures,  heat  transfer  in  mass  flow  sensors,  and 
unsteady  pressure  fields  around  microvalves,  for  example,  must 
be  understood  if  these  devices  are  to  be  effectively  designed 
and  optimized. 

Unfortunately,  Navier-Stokes-based  computational  fluid  dy¬ 
namics  (CFD)  techniques  are  often  inaccurate  when  applied  to 
MEMS.  This  inaccuracy  stems  from  their  calculation  of  molecu¬ 
lar  transport  effects,  such  as  viscous  dissipation  and  thermal 
conduction,  from  bulk  flow  quantities,  such  as  mean  velocity 
and  temperature.  This  approximation  of  micro-scale  phenomena 
with  macro-scale  information  fails  as  the  characteristic  length 
of  the  flow  gradients  {£)  approaches  the  average  distance  trav¬ 
eled  by  molecules  between  collisions  (the  mean  free  path,  X.). 
The  ratio  of  these  quantities  is  known  as  the  Knudsen  number 
=  X/i?)  and  is  used  to  indicate  the  degree  of  flow  rarefac¬ 
tion.  For  Kn  <  0.01,  the  flow  is  considered  to  be  in  the  ‘contin¬ 
uum’  regime  and  the  Navier-Stokes  equations  are  applicable  in 
their  common  form.  As  Kn  increases,  the  flow  moves  through 
the  “slip-flow”  (0.01  <  Kn  <  0.1)  and  “transition”  (0.1  < 
Kn  <  3)  regimes  and  finally  enters  the  “free-molecular”  (Kn 
>3)  regime,  each  suggesting  a  particular  type  of  analysis 
(Schaff  and  Chambre,  1958). 

The  Knudsen  number  of  many  MEMS-related  flows  is  driven 
from  the  continuum  regime  by  an  extremely  small  feature  size, 
which  is  often  comparable  to  \,  even  for  air  at  standard  condi¬ 
tions.  Noncontinuum  effects,  neglected  in  traditional  analyses. 
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may  therefore  significantly  affect  device  performance.  As  a  re¬ 
sult,  new  models  and  analysis  techniques  must  be  developed  to 
correctly  predict  the  fluid  behavior  in  and  around  MEMS. 

For  numerical  analyses,  the  direct  simulation  Monte  Carlo 
(DSMC)  method  (Bird,  1994)  offers  an  alternative  to  tradi¬ 
tional  CFD  that  retains  its  validity  at  high  Knudsen  numbers. 
Commonly  applied  to  reentry  vehicles,  this  technique  makes  no 
continuum  assumption.  Instead,  it  models  the  flow  as  it  physi¬ 
cally  exists:  a  collection  of  discrete  particles,  each  with  a  posi¬ 
tion,  a  velocity,  an  internal  energy,  a  species  identity,  etc.  These 
particles  are  moved  and  allowed  to  interact  with  the  domain 
boundaries  in  small  time  steps  during  the  calculation.  Intermo- 
lecular  collisions  are  all  performed  on  a  probabilistic  basis  at 
the  end  of  each  time  step  to  minimize  computational  work. 
Macroscopic  quantities,  such  as  flow  speed  and  temperature, 
are  then  obtained  by  sampling  the  microscopic  state  of  all  parti¬ 
cles  in  the  region  of  interest. 

Algorithm 

The  DSMC  code  developed  for  this  investigation  is  written 
for  unstructured  grids  to  provide  flexibility  for  complex  geome¬ 
tries.  Particle  information  is  stored  locally  with  each  cell,  fol¬ 
lowing  Dietrich  and  Boyd  (1994),  to  take  advantage  of  the 
cache  hardware  of  workstation-class  computers  and  to  prepare 
for  a  planned  transition  to  parallel  machines.  Bird’s  No  Time 
Counter  (NTC)  scheme  ( 1989)  governs  collision  pair  selection, 
with  collision  cross-sections  calculated  according  to  his  Vari¬ 
able  Hard  Sphere  (VHS)  model  (1981). 

Particle  movement  and  current-cell  identification  are  per¬ 
formed  simultaneously  with  a  trajectory- tracing  scheme  similar 
to  that  proposed  by  Dietrich  ( 1991 ).  In  this  scheme,  a  particle 
is  moved  along  its  trajectory  until  it  contacts  a  face  of  its  current 
cell.  It  is  then  passed  to  another  cell,  reflected  from  a  solid 
surface,  or  allowed  to  leave  the  calculation,  as  appropriate. 

Particles  impinging  on  solid  walls  are  reflected  diffusely  with 
full  thermal  and  momentum  accommodation.  Those  encounter¬ 
ing  inflow /outflow  faces  are  simply  allowed  to  leave  the  do¬ 
main.  The  cell-based  storage  of  particle  information  simplifies 
the  latter  operation;  impinging  particles  are  simply  transferred 
out  of  their  current  cells  without  passing  them  to  adjacent  cells. 

At  I/O  faces,  particles  must  be  introduced  to  represent  the 
influence  of  fluid  outside  the  domain.  In  accord  with  standard 
CFD  practice,  their  quantity  and  velocity  distribution  is  deter¬ 
mined  by  considering  the  flow  “characteristics”  (cf.  Hirsch, 
1990),  or  lines  along  which  certain  variables  remain  constant; 
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namely  the  entropy,  the  transverse  speed,  and  the  Riemann 
invariants,  which  are  given  (-for  isentropic  flow)  by: 

T  2a  ,  _ 

J±  =  u± - -  ,  (1) 

7  -  1 

where  u  is  the  streamwise  speed,  a  is  the  speed  of  sound,  and 
7  is  the  ratio  of  specific  heats. 

The  paths  of  these  characteristics  dictate  how  much  informa¬ 
tion  may  be  specified  at  the  boundary  a  priori,  with  the  remain¬ 
der  obtained  from  the  instantaneous  state  inside  the  domain.  In 
a  supersonic  flow,  all  characteristics  point  downstream.  In  a 
subsonic  flow,  the  J-  characteristic  points  upstream  and  all 
others  point  downstream. 

Most  DSMC  applications  to  date  have  involved  high  speed 
flows.  This  greatly  simplifies  I/O  boundary  treatment  because 
all  variables  at  the  inlet  may  be  specified  by  the  user.  At  the 
outlet,  where  it  would  be  necessary  to  determine  all  variables 
from  the  flowfield,  no  treatment  is  required  because  an  insig¬ 
nificant  number  of  particles  travel  upstream  into  the  domain, 
making  it  possible  to  neglect  the  particle  introduction  step  alto¬ 
gether. 

Most  MEMS,  however,  involve  low  speed  flows.  It  is  there¬ 
fore  necessary  to  determine  some  variables  at  both  inlet  and 
outlet  faces  from  inside  the  flowfield.  All  cases  discussed  in 
this  article  are  intended  to  represent  portions  of  long  devices  far 
from  their  pressure  reservoirs.  To  accomplish  this,  the  pressure, 
temperature,  and  transverse  speed  are  specified  at  itic  inlet, 
leaving  the  streamwise  velocity  to  be  obtained  from  inside  the 
domain.  Only  the  pressure  is  specified  at  the  outlet  and  the 
remaining  variables  are  determined  from  inside  the  domain. 

In  a  continuum  code,  this  operation  is  straightforward;  the 
node  just  inside  the  boundary  is  simply  queried  for  the  necessary 
information.  Unfortunately,  DSMC’s  statistical  nature  compli¬ 
cates  this  step  because  there  may  be  as  few  as  20  particles  in 
a  given  cell  at  a  given  time.  Determining  its  macroscopic  state 
from  an  instantaneous  sample  therefore  yields  unacceptable  sta¬ 
tistical  scatter.  Two  immediate  options  exist  for  a  steady  flow: 
neighboring  cells  may  be  included  in  the  instantaneous  sample 
or  it  may  be  replaced  with  a  time  average.  The  proper  applica¬ 
tion  of  the  former  option  requires  identifying  neighboring  cells 
with  states  “close-enough”  to  the  cell  in  question  to  yield  a 
meaningful  spatial  average.  The  latter  option  involves  choosing 
a  time  averaging  method  that  yields  a  sufficiently  accurate  esti¬ 
mate  but  still  allows  a  boundary  to  be  reasonably  responsive  to 
changes  in  the  flow  as  it  seeks  its  steady  state. 

The  latter  option  is  implemented  in  the  current  code.  An 
inflow /outflow  cell’s  particles  are  sampled  after  a  movement 
step  is  completed,  incoming  particles  are  introduced  at  bound¬ 
aries,  and  collisions  are  performed.  A  weighted  average  is  then 
taken  between  this  result  and  a  running  value  collected  from 
previous  time  steps.  The  weighting  of  the  instantaneous  state 
may  be  varied  to  balance  the  speed  of  convergence  with  the 
accuracy  of  the  running  estimate.  A  weight  of  ^  was  chosen 
for  the  cases  presented  in  this  article. 

Results 

Slip  Flow  Regime  Micro-Channel.  The  first  case  explored 
with  this  code  was  a  steady,  low-speed,  2-D  flow  through  a 
micro-channel  with  an  outlet  Knudsen  number  of  0.05,  placing 
it  in  the  slip  flow  regime.  This  geometry  is  similar  to  those 
investigated  experimentally  by  Harley  et  al.  (1995),  Pong  et 
al.  (1994),  and  Arkilic  et  al.  (1994)  and  numerically  (with 
spectral  element  methods)  by  Beskok  and  Kamiadakis  ( 1994) . 
This  is,  historically,  an  important  canonical  case  for  determining 
the  effect  of  rarefaction  on  the  transport  terms  in  the  Navier- 
Stokes  equations.  It  is  also  representative  of  flows  along  the 
narrow  passages  often  found  in  MEMS  devices,  such  as  the 
space  under  the  floating  plate  of  a  shear  stress  sensor  or  acceler- 
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Fig.  1  Comparison  of  computed  and  analytical  pressure  distributions 
for  a  micro-channel  in  the  slip  flow  regime.  Excellent  agreement  is  ob¬ 
tained  between  code  and  theory  including  the  nonlinear  pressure  distri¬ 
bution  resulting  from  compressibility. 


ometer,  which  is  typically  only  one  micron  high  but  hundreds 
of  microns  in  breadth  and  depth  (Padmanabhan  et  al.,  1995). 
In  addition,  it  serves  as  a  convenient  calibration  case  to  assess 
the  accuracy  of  the  numerical  algorithm  because  analytical  solu¬ 
tions  have  been  developed  for  this  geometry. 

In  one  such  effort,  Arkilic  et  al.  show  that  the  Navier-Stokes 
equations  may  be  solved  analytically  for  a  long,  high  aspect- 
ratio,  isothermal  channel  in  the  slip  flow  regime  if  the  boundary 
conditions  are  modified  to  include  a  Kn-dependent  streamwise 
velocity  (slip)  at  the  wall,  given  (for  diffuse  reflection)  by: 

=  (2) 

dy  ^-all 

where  u  is  the  streamwise  velocity,  Kn  is  the  local  Knudsen 
number  and  y  is  the  transverse  coordinate,  which  has  its  zero 
at  the  channel  centerline. 

Through  this  analysis,  an  expression  may  be  obtained  for  the 
pressure  distribution  as  a  function  of  streamwise  channel  loca¬ 
tion  and  overall  pressure  ratio: 

f(x)  =  -6  Kn, 

+  y(6  Kn„  +  -  £  [iff  -  1)  +  12  Kn,if  -  1)]  (3) 

where  f(x)  and  %  are  the  local  and  inlet  pressures,  respectively, 
normalized  by  the  outlet  value,  Kn,  is  the  outlet  Knudsen  num¬ 
ber,  X  is  the  streamwise  coordinate,  and  L  is  the  channel  length. 

The  distribution  predicted  by  Eq.  (3)  may  be  compared  to  a 
DSMC  result  as  a  test  of  both  theory  and  code.  Such  a  compari¬ 
son  is  presented  in  Fig.  1  for  a  31.14  X  1.04  fim  channel  filled 
with  nitrogen  at  a  pressure  ratio  of  2.47.  A  good  agreement  is 
obtained  (max  error  =  1.5  percent),  including  the  pressure 
curve  non-linearity.  This  nonlinearity  is  a  result  of  the  large 
down-channel  pressure  drop,  which  causes  a  significant  density 
variation  (the  flow  is  essentially  isothermal).  This  was  observed 
experimentally  by  Pong  et  al.  in  silicon  microchannels  with 
integral  pressure  taps  and  the  nonlinearity  was  found  to  become 
more  pronounced  as  the  pressure  ratio  increased.  Interestingly, 
the  form  of  disagreement  between  frie_analytic  results  and 
DSMC  presented  here  mirrors  that  found  by  Beskok  et  al. 
( 1996)  between  their  spectral  code  (with  modified  wall  bound¬ 
ary  conditions)  and  a  DSMC  code  of  Bird  (1994). 
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Fig.  2  Theoretical  (— )  and  computed  ( — )  streamwise  velocity  distribu¬ 
tions  (in  m/s)  fora  micro-channel  in  the  slip-flow  regime.  A  good  agree¬ 
ment  is  obtained,  including  the  flow  acceleration  due  to  compressibility 
and  the  non-zero  velocity  at  the  wall  due  to  rarefaction. 


A  theoretical  expression  for  the  streamwise  velocity  distribu¬ 
tion  was  also  developed  by  Arkilic  et  al.: 


dx 


-  H^Kn 


(4) 


where  ^  is  the  coefficient  of  viscosity,  p  is  the  pressure,  and  H 
is  the  channel  height. 

This  equation  is  plotted  in  Fig,  2,  along  with  the  DSMC 
result,  for  the  geometry  and  conditions  used  above.  An  excellent 
agreement  between  theory  and  code  is  apparent,  despite  the 
statistical  fluctuation  present  in  the  DSMC  contours.  Unfortu¬ 
nately,  this  fluctuation  is  an  inherent  feature  of  low-speed 
DSMC  output  because  the  macroscopic  flow  speed  is  small 
compared  to  the  mean  particle  thermal  speed.  Calculating  the 
flow  speed  is  therefore  a  matter  of  isolating  a  small-amplitude 
signal  superimposed  on  large-amplitude  noise;  an  operation 
which  requires  an  enormous  number  of  samples. 

Several  interesting  features  common  in  micro-channel  flows 
are  visible  in  this  figure.  First,  the  fluid  accelerates  as  it  moves 
down  the  channel,  unlike  in  the  familiar  Poiseuille  result.  This 
is  another  consequence  of  the  density  drop  caused  by  the  de¬ 
creasing  pressure  in  the  streamwise  direction;  the  mean  stream- 
wise  velocity  must  increase  to  maintain  a  constant  mass  flow. 
Second,  the  velocity  at  the  walls  is  nonzero  and  increases  with 
increasing  x-coordinate.  This  is  the  aforementioned  “slip  flow,” 
which,  by  Eq.  (2),  is  negligible  for  continuum  flows  due  to 
their  very  small  Knudsen  numbers.  The  increase  in  slip  velocity 
down  the  channel  is  caused  by  growth  in  both  Kn  (from  the 
decreasing  pressure)  and  velocity  gradient  at  the  wall  (from  the 
accelerating  flow).  This  is  contrary  to  the  behavior  very  close 
to  the  entrance  observed  by  Beskok  and  Kamiadakis  (1994) 
and  Oh  et  al.  (1995),  which  is  absent  from  this  work  because 
the  boundary  conditions  were  formulated  to  eliminate  entrance 
effects. 

A  further  comparison  to  the  theoretical  analysis  of  Arkilic  et 
al.  may  be  obtained  by  normalizing  the  velocity  distribution  of 
Eq.  (4)  by  the  average  velocity  at  a  given  x-location,  Uave, 
obtained  by  integrating  u  from  the  lower  wall  to  the  upper  wall 
and  dividing  by  H.  After  some  rearrangement,  this  yields: 


~Kn  + 


u  _  I 
^ave  4 


(5) 


The  left  side  of  this  equation,  which  will  be  referred  to  as  the 


“similarity  speed,”  is  a  function  of  x  and  y,  while  the  right 
is  a  function  of  y  alone.  Consequently,  if  the  slip-flow  analysis 
holds,  calculating  the  similarity  speed  using  the  local  Kn(x) 
and  m(x,  y)  will  eliminate  the  x-dependence  of  the  velocity 
distribution  shown  in  Fig,  2,  collapsing  the  profile  at  each  x- 
location  to  a  single  parabola  given  by  the  right  side  of  Eq.  (5). 

This  assertion  was  tested  by  computing  a  similarity  speed 
distribution  from  the  DSMC  output.  As  predicted  by  the  analy¬ 
sis,  it  was  found  that  the  down-channel  variation  of  the  stream- 
wise  velocity  profile  seen  in  Fig.  2  collapsed  to  a  constant 
similarity  speed  profile.  It  was  therefore  concluded  that  the 
similarity  assertion  expressed  in  Eq.  (5)  indeed  holds  for  the 
slip-flow  channel. 

Overall,  excellent  agreement  was  obtained  between  the  ana¬ 
lytical  solution  of  Arkilic  et  al.  and  the  DSMC  results.  This 
supports  the  accuracy  of  both  techniques.  For  the  DSMC  code, 
however,  this  a  just  a  convenient  validation  case;  many  more 
interesting  flows,  for  which  there  are  no  reliable  analytical  solu¬ 
tions,  may  be  easily  treated  with  this  method.  The  remaining 
cases  presented  in  this  article  are  intended  to  demonstrate  this 
capability. 

Transition  Regime  Micro-Channel.  One  of  DSMC’s 
greatest  strengths  is  its  validity  for  dilute  gases  in  all  Knudsen 
number  regimes.  By  far,  the  most  unexplored  of  these  is  the 
transition  regime,  0.1  <  Kn  <  3.  Obtaining  solutions  at  these 
Knudsen  numbers  is  very  difficult  because  the  approximation 
of  transport  terms  based  on  macroscopic  quantities  becomes 
unacceptably  inaccurate,  precluding  the  use  of  the  Navier- 
Stokes  equations,  even  with  the  slip-flow  boundary  condition 
of  Eq.  (2).  Collisions  are  still  important,  however,  so  the  colli- 
sionless  Boltzmann  equation  is  not  yet  an  option. 

Because  DSMC  is  a  direct  method,  rather  than  a  numerical 
solution  of  a  discretized  system  of  equations,  it  is  a  veiy  attrac¬ 
tive  tool  for  investigating  the  transition  regime.  This  direct  na¬ 
ture  also  makes  it  one  of  the  few  tractable  techniques  which  is 
uniformly  valid  in  mixed  Kn-regime  flows.  These  are  very  use¬ 
ful  features  because,  due  to  the  aforementioned  problems,  rela¬ 
tively  little  is  known  about  these  cases.  Such  knowledge  is 
critical  to  MEMS  designers,  however,  because  many  devices 
contain  flows  characterized  by  these  difficulties. 

To  highlight  the  changes  that  occur  as  Kn  enters  the  transition 
regime,  a  micro-channel  is  again  used  as  the  sample  case.  This 
time,  however,  the  analytical  solution  of  Arkilic  et  al.  serves 
only  as  a  convenient  reference,  allowing  transition  regime  re¬ 
sults  to  be  compared  to  continuum  (Kn  =  0.0)  and  slip-flow 
predictions;  significant  disagreement  is  expected.  The  channel 
treated  was  6.61  X  0.11  jjm,  with  a  pressure  ratio  of  4.2  and 
an  outlet  Knudsen  number  of  0.44.  The  working  fluid  was, 
again,  nitrogen. 

Proceeding  along  the  same  path  of  comparison  used  in  the 
previous  case,  the  computed  pressure  distribution  is  plotted  with 
the  analytical  prediction  in  Fig.  3.  The  continuum  curve  (Kn 
=  0.0)  is  also  shown  for  reference.  As  expected,  the  excellent 
agreement  between  DSMC  and  theory  obtained  for  the  slip- 
flow  case  (Fig.  1)  is  no  longer  present.  The  error  between  the 
curves  has  grown  from  less  than  2  to  almost  5  percent.  The 
form  of  this  disagreement  is  also  significant:  the  computed  curve 
is  more  linear  than  its  analytical  counterpart.  A  trend  of  increas¬ 
ing  pressure  curve  linearity  with  increasing  rarefaction  is  there¬ 
fore  established  by  the  relative  shape  of  the  continuum,  slip- 
flow,  and  transition  curves. 

This  trend  is  contr^  to  the  experimental  observations  of 
Pong  et  al.  (1994)  (Fig.  7),  which  were  obtained  by  running 
a  channel  with  nitrogen  and  then  with  helium  at  the  same  pres¬ 
sure  ratio.  It  agrees,  however,  with  the  conclusions  of  Beskok 
et  al.  ( 1995) .  The  source  of  this  discrepancy  between  analytical, 
DSMC,  and  spectral  results  with  this  experimental  data  has  yet 
to  be  identified. 
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Fig.  3  DSMC-compifted  pressure  distribution  for  a  transition  regime 
micro-channel  compared  to  analytical  (slip-flow)  solutions  for  the  same 
Kn  and  for  continuum  flow  (Kn  -  0.0).  A  trend  of  increasing  pressure 
curve  linearity  with  increasing  rarefaction  is  apparent. 


Fig.  5  Computed  streamwise  velocity  distribution  for  the  channel  of 
Figure  4.  The  magnitude  of  both  the  slip  and  the  local  maximum  velocity 
are  greater  than  the  analytical  prediction  by  as  much  as  40%. 


The  streamwise  velocity  distribution  predicted  by  Eq.  (4)  for 
this  channel  is  presented  in  Fig.  4.  Comparing  this  to  Fig.  2,  it 
is  evident  that  Eq.  (4)  yields  a  much  flatter  profile  for  the  higher 
Kn  case.  This  may  be  attributed  to  the  last  term  in  Eq.  (4), 
which  is  constant  across  the  channel  and  proportional  to  Kn. 
The  Knudsen  number  is  an  order  of  magnitude  larger  in  this 
case,  so  the  constant  term  has  a  much  stronger  influence  on  the 
shape  of  the  distribution. 

Upon  plotting  the  DSMC  result  for  comparison  in  Fig.  5,  it 
again  becomes  evident  that  the  assumptions  supporting  Eq.  (4) 
are  beginning  to  fail.  Both  the  slip  flow  and  maximum  speeds 
at  a  given  x-location  are  higher  than  predicted  by  as  much  as 
40  percent.  This  allows  the  channel  to  support  a  much  larger 
mass  flow  than  predicted  by  the  slip-flow  solution,  which  is,  in 
turn,  a  larger  mass  flow  than  predicted  by  the  no-slip  solution. 

Using  the  similarity  expression  developed  in  the  previous 
section,  the  DSMC  result  may  be  used  to  determine  the  Kn 
limit  for  the  slip-flow  solution,  commonly  given  as  Kn  <  0.1. 
Because  each  x-location  has  both  a  similarity  profile  and  a 
Knudsen  number  associated  with  it,  the  Kn  at  which  the  slip- 
flow  analysis  fails  may  be  determined  by  finding  the  point  where 
the  DSMC  and  analytical  profiles  begin  to  disagree.  Toward 
this  end,  Fig.  6  contains  the  computed  maximum  and  minimum 


similarity  speeds  at  each  x-location,  plotted  with  the  analytical 
predictions  (which  form,  by  definition,  straight  lines).  The  Kn 
distribution  was  then  overlaid  to  facilitate  determining  its  value 
when  the  slip-flow  analysis  fails. 

It  is  evident  from  this  figure  that  the  slip  flow  analysis  begins 
to  fail  at  approximately  Kn  =  0.15.  This  limit  may  be  under¬ 
stood  if  the  boundary  condition  of  Eq.  (2)  is  viewed  as  an 
expansion  of  the  wall  velocity  in  powers  of  Kn.  The  no-slip 
condition  and  Eq.  (2)  are  then  the  zeroth  and  first  order  approxi¬ 
mations,  respectively.  It  is  therefore  logical  that  the  neglected 
second  and  higher-order  terms  would  begin  to  significantly  af¬ 
fect  the  result  when  Kn  exceeds  approximately  0.1.  A  second- 
order  accurate  boundary  condition  in  terms  of  the  continuum 
variables  is  presented  in  Beskok  et  al.  (1994).  It  should  be 
remembered,  however,  that  the  Navier-Stokes  equations  them¬ 
selves  are  only  strictly  valid  to  first  order  in  ICn. 

Supersonic  Micro-Nozzle.  The  final  case  presented  in  this 
article  is  a  supersonic  micro-nozzle.  In  relation  to  the  previous 
cases,  this  geometry  may  be  viewed  as  a  channel  whose  upper 
and  lower  walls  form  a  parabolic  contraction /expansion.  The 
expansion  region  of  similar  nozzles  was  treated  with  DSMC  by 
Zelesnik  et  al.  (1994)  with  application  to  satellite  station-keep¬ 
ing  rocket  motor  design.  The  nozzle  treated  here  is  much  smaller 
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Fig,  4  Theoretical  streamwise  velocity  distribution  for  a  transition-re¬ 
gime  micro  channel  calculated  using  the  slip-flow  analysis 
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Fig.  6  DSMC  (— )  and  analytical  {--)  similarity  speeds  with  a  Knudsen 
number  (-*-)  overlay.  The  point  at  which  the  computed  curve  begins  to 
diverge  from  the  prediction  corresponds  well  with  the  commonly-stated 
Kn  for  the  onset  of  the  transition  regime  (0.1 ) . 
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however — only  15.4  pbvn  high  at  the  throat  and  92.6  long.  It 
is  also  2-D,  rather  than  axisymmetric,  which  is  a  more  attractive 
configuration  for  a  MEMS  device  due  to  the  inherently  planar 
nature  of  photolithographic  manufacturing  techniques. 

Such  a  nozzle  may  be  considered  quasi- ID  if  its  area  change 
is  sufficiently  gradual.  The  analytical  solution  of  Arkilic  et  al. 
could  then  be  used,  with  appropriate  modifications  for  the  x- 
dependent  area.  In  the  case  discussed  here,  however,  the  total 
length  is  only  six  times  the  throat  height.  In  addition,  the  large 
density  variations  brought  about  by  the  rapid  expansion  may 
cause  the  Kn-regime  to  change  at  one  or  more  streamwise  posi¬ 
tions. 

These  factors  make  both  analytical  and  continuum-based  nu¬ 
merical  treatment  of  this  geometry  prohibitively  difficult.  None¬ 
theless,  nozzles  such  as  these  may  play  important  roles  in 
MEMS  devices  like  micro-rocket  thrusters  and  micro-gas  tur¬ 
bine  generators,  for  example.  Investigating  their  behavior  is 
therefore  a  valuable  task  for  which  DSMC  is  well-suited. 

Due  to  the  large  pressure  variation  across  the  nozzle  and 
DSMC’s  requirement  that  a  cell’s  linear  dimensions  be  propor¬ 
tional  to  the  local  mean  free  path,  26,758  cells  were  required 
to  properly  discretize  the  domain.  The  inlet  was  charged  with 
helium  at  1  atmosphere  and  a  ‘vacuum’  outlet  condition  was 
enforced.  This  condition  was  implemented  by  simply  removing 
from  the  simulation  any  particle  crossing  the  outlet  boundary 
and  refraining  from  introducing  new  particles  at  these  faces. 
The  resulting  pressure  at  the  outlet  plane  was  approximately 
0.04  atmospheres,  yielding  an  outlet  Knudsen  number,  based 
on  passage  height,  of  0.03  and  a  pressure  ratio  of  24.  Approxi¬ 
mately  550,000  particles  were  present  on  the  grid  at  steady 
state. 

The  Mach  number  distribution  for  this  nozzle  is  shown  in 
Fig.  7.  Qualitatively,  the  contours  agree  quite  well  with  the 
results  of  Zelesnik  et  al.  In  addition,  a  number  of  interesting 
features  are  visible.  First,  contrary  to  the  inviscid,  quasi- ID 
result,  the  sonic  point  is  actually  slightly  downstream  of  the 
throat.  This  is  possible  because  the  viscous  effects  are  strong 
enough  to  overcome  the  deceleration  due  to  the  diverging  shape. 
Second,  a  Mach  number  of  2.4  is  reached;  a  considerable  speed 
for  such  a  short,  narrow  nozzle.  Finally,  the  slip  flow  speed  is 
substantial,  exceeding  Mach  0,5  near  the  outlet. 

The  temperature  distribution  for  this  micro-nozzle  is  shown 
in  Fig.  8.  Again,  a  good  qualitative  agreement  with  the  results 
of  Zelesnik  et  al.  is  obtained.  Clearly,  unlike  in  the  previous 
two  cases,  this  flow  cannot  be  considered  isothermal.  A  strong 
temperature  gradient  exists  in  both  the  streamwise  and  trans¬ 
verse  directions,  creating  yet  another  obstacle  to  analytical  treat¬ 
ment.  The  presence  of  such  strong  temperature  variation  is  a 
striking  feature  when  the  diminutive  dimensions  of  the  nozzle 
are  considered.  Though  the  walls  are  isothermal  with  full  energy 
accomodation  and  only  approximately  30  local  mean  free  paths 
apart  at  the  exit,  the  fluid  is  still  able  to  realize  a  substantial 
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Fig.  7  Mach  number  distribution  for  a  supersonic  micro-nozzle.  A  maxi¬ 
mum  Mach  number  of  2.4  is  obtained  from  a  pressure  ratio  of  24.  The 
slip  speed  exceeds  Mach  0.5  near  the  outlet. 


Fig.  8  Temperature  contours  for  a  supersonic  micro-nozzle,  normalized 
by  the  inlet  value.  Significant  temperature  gradients  are  present,  despite 
the  fact  that  the  walls  are  isothermal  (at  the  inlet  temp)  with  full  thermal 
accommodation. 


reduction  in  temperature.  The  effect  of  rarefaction  has  therefore 
been  to  significantly  reduce  the  thermal  communication  between 
the  wall  and  the  fluid.  This  assertion  is  further  supported  by 
noting  the  large  thermal  slip  at  the  wall. 

Conclusion 

The  direct  simulation  Monte  Carlo  method  has  proven  to  be 
a  valuable  tool  for  investigating  the  behavior  of  flows  which 
are  considered  ‘rarefied’  due  to  miniaturization.  This  article  has 
presented  only  a  small  subset  of  the  many  cases  of  interest  to 
both  the  fluid-dynamicist  and  the  MEMS  designer  which  are 
easily  treated  with  this  method. 

The  unique  features  of  rarefied  gas  flows  can  significantly 
affect  MEMS  operation.  The  larger  mass  flow  rate  for  a  given 
geometry  and  inlet  conditions  must  be  considered  when  de¬ 
signing  control  systems  for  micro-chemical  reactors  and  the 
decreased  thermal  cormnunication  between  the  fluid  and  its 
surroundings  must  be  considered  for  applications  involving 
temperature  measuring  devices  or  heaters,  for  example.  Unfor¬ 
tunately,  these  effects  are  neglected  in  many  common  design 
tools. 

DSMC,  therefore,  has  great  promise  for  improving  the  design 
and  optimization  of  MEMS.  Its  ability  to  calculate  in  any  of 
the  four  Knudsen  number  regimes  without  modification  is  an 
important  feature  for  this  task.  This  is  especially  valuable  for 
flows  with  regions  in  different  regimes.  The  addition  of  an 
unstructured  grid  capability  and  a  trajectory-tracing  movement 
scheme  enables  the  code  to  handle  arbitrary  geometries.  This 
is  a  valuable  asset  when  analyzing  the  complex  structures  found 
in  many  devices.  In  addition,  it  is  quite  straightforward  to  in¬ 
clude  flow  features  such  as  chemical  reactions,  multi-species 
mixing  and  particle  transport  in  the  code  due  to  its  modular 
nature.  Finally,  its  cell-based  structure  makes  it  well-suited  for 
parallel  computation,  which  is  an  increasingly-important  attri¬ 
bute  for  a  large-scale  numerical  scheme. 
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A  DSMC  investigation  of  flows  related  to  micro-electro-mechanical  systems  (MEMS)  is  detailed.  This  effort  is 
aimed  at  increasing  our  understanding  of  rarefled  gas  behavior  to  facilitate  the  design  and  optimization  of  micro¬ 
devices.  The  code  written  for  this  work  employs  an  unstructured  grid  and  a  trajectory-tracing  particle  movement 
scheme.  Its  results  for  slip-flow  and  transition  regime  micro-channels  and  a  micro-nozzle  are  presented.  The  slip  flow 
micro-channel  output  Is  compared  to  an  analytical  solution  of  the  Navier-Stokes  equations  with  a  first-order  (in  Kn) 
slip  boundary  condition  at  the  wail.  Excellent  agreement  is  found.  A  transition  regime  micro-channel  is  also 
compared  to  the  analytical  solution  in  order  to  investigate  the  effects  of  further  rarefaction  on  flow  behavior.  A 
higher  mass  flow  rate  and  a  lower  degree  of  pressure  distribution  nonlinearity  is  observed.  The  micro-nozzle 
produces  a  maximum  Mach  number  of  2.4  from  a  pressure  ratio  of  24.  Though  its  walls  are  held  at  the  reference 
temperature  and  full  thermal  accommodation  is  implemented,  a  large  temperature  gradient  is  established  in  both  the 
streamwise  and  transverse  directions.  This  allows  the  flow  to  reach  a  minimum  temperature  which  is  less  than  40% 
of  the  wall  value. 


L  Introduction 

MICRO-electro-mechanical  systems  (MEMS)  are 
the  subject  of  increasingly  active  research  in  a 
widening  field  of  disciplines.  These  are  very  small 
(micron-scale)  sensors  and  actuators  manufactured  with 
techniques  similar  to  those  used  for  integrated  circuit 
(IC)  chips.  Applications  for  these  devices  range  from 
consumer  products  (air-bag  triggers,  micro-mirror  dis¬ 
plays),  to  industrial  and  medical  tools  (micro- valves, 
micro- motors),  to  instrumentation  (micro  pressure  sen¬ 
sors,  micro  shear-stress  sensors).  ^ 

MEMS  have  many  advantages  over  their  macro¬ 
scale  counterparts,  where  such  counterparts  even  exist. 
First,  because  these  devices  are  fabricated  in  a  manner 
similar  to  IC  chips,  they  are  extremely  inexpensive  to 
manufacture  in  large  quantities.  Second,  the  technology 
for  such  production  is  quite  mature.  Very  precise  specifi¬ 
cation  of  the  geometry,  far  beyond  that  possible  with 
macro-scale  fabrication  techniques,  is  therefore  routine 
and  a  high  degree  of  control  over  material  properties  is 
possible.  Third,  their  small  size  and  mass  make  them 
attractive  where  space  is  at  a  premium  or  weight  is  lim¬ 
ited.  Finally,  their  minimal  inertia  allows  them  to  react 
very  quickly,  enabling  the  creation  of  extremely  fast 
actuators  and  sensors  with  previously  unthinkable  fre¬ 
quency  response. 
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The  small  size  of  MEMS  poses  unique  challenges  in 
the  design  phase,  however.  While  the  mechanical  prop¬ 
erties  of  micromachined  materials  are  reasonably  well- 
studied,  fluid  effects  at  micron  scales  are  not.  These 
effects,  such  as  film  damping  of  resonant  structures,  heat 
transfer  in  mass  flow  sensors,  and  unsteady  pressure 
fields  around  microvalves,  for  example,  must  be  under¬ 
stood  if  the  full  potential  of  these  devices  is  to  be  real¬ 
ized. 

Unfortunately,  traditional  CFD  techniques  are  often 
inaccurate  for  analyzing  MEMS.  At  fault  are  the  trans¬ 
port  terms  in  the  Navier-Stokes  equations  which  are  cal¬ 
culated  from  macro-scale  flow  quantities,  such  as 
temperature  and  mean  velocity.  This  approximation  of 
micro-scale  phenomena  with  macro-scale  information 
fails  as  the  characteristic  length  of  the  flow  gradients  (L) 
approaches  the  average  distance  traveled  by  molecules 
between  collisions  (the  mean  free  path,  X).  The  ratio  of 
these  quantities  is  known  as  the  Knudsen  number  (Kn  = 
X/  L)  and  is  used  to  indicate  the  degree  of  flow  rarefac¬ 
tion.  The  Navier-Stokes  equations  ignore  rarefaction 
effects  and  are  therefore  only  strictly  accurate  at  a  van¬ 
ishingly  small  Kn. 

As  is  the  case  for  other  nondimensional  numbers  in 
fluid  dynamics,  such  as  the  Mach  and  Reynolds  num¬ 
bers,  the  type  of  analysis  appropriate  for  a  particular 
flow  is  dictated  by  its  Knudsen  number.  The  Kn  domain 
(0  <  Kn  <  oo)  is  often  divided  into  four  flow  regimes.- 
For  Kn  <0.01,  known  as  the  ‘continuum’  regime,  the 
Navier-Stokes  equations,  as  commonly  expressed,  are 
applied.  For  0.01  <  Kn  <  0.1,  known  as  the  ‘slip-flow’ 
regime,  the  Navier-Stokes  equations  are  applied  with 
the  commonly-used,  no-slip  wall  boundary  condition 
replaced  with  a  slip-flow  condition  (to  be  discussed 
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later).  For  0.1<Kn<3»  known  as  the  ‘transition’ 
regime,  the  flow  is  too  rarefied  for  Navier-Stokes-based 
analysis,  but  not  rarefied  enough  to  apply  the  collision¬ 
less  Boltzmann  equation.  The  full  Boltzmann  equation 
is  therefore  prescribed.  For  Kn  >  3,  known  as  the  ‘free 
molecular’  regime,  the  flow  is  sufficiently  rarefied  to 
allow  molecular  collisions  to  be  completely  neglected  in 
analyses.  The  collisionless  Boltzmann  equation  is  there¬ 
fore  applied. 

For  many  MEMS,  Kn  is  driven  from  the  continuum 
regime  by  an  extremely  small  feature  size,  which  is 
often  comparable  to  K  even  at  standard  conditions.  The 
gap  between  the  sensing  plate  and  the  substrate  on  a 
floating  element  shear  stress  sensor,  for  example,  is  typ¬ 
ically  1-2  microns.  The  mean  free  path  of  air  at  standard 
conditions  is  approximately  60nm.  This  places  Kn  in  the 
slip  flow  regime,  or  even  the  transition  regime  for 
lighter  gases  or  different  conditions.  Non-continuum 
effects  in  the  gap,  neglected  in  traditional  analyses,  may 
therefore  have  a  significant  impact  on  sensor  operation. 
As  a  result,  new  models  and  techniques  must  be  devel¬ 
oped  to  correctly  describe  the  behavior  of  the  fluid  in 
and  around  these  devices. 

For  numerical  analyses,  the  Direct  Simulation 
Monte  Carlo  (DSMC)  method  offers  an  alternative  to 
traditional  CFD  algorithms  which  retains  its  validity  at 
high  Knudsen  numbers.  Commonly  applied  to  space 
vehicles,  this  technique  makes  no  continuum  assump¬ 
tion;  it  models  the  flow  as  it  physically  exists:  a  collec¬ 
tion  of  discrete  particles,  each  with  a  position,  a 
velocity,  and,  if  required,  an  internal  energy.  These  par¬ 
ticles  are  moved  and  allowed  to  interact  with  each  other 
and  the  domain  boundaries  as  appropriate  during  the 
calculation.  Macroscopic  quantities,  such  as  flow  speed 
and  temperature,  are  then  obtained  by  sampling  the 
microscopic  state  of  all  particles  in  the  region  of  inter¬ 
est. 

11.  Algorithm 

MEMS  devices  typically  contain  a  variety  of  com¬ 
plex  geometries.  To  provide  the  flexibility  necessary  to 
treat  these  cases,  the  DSMC  code  developed  for  this 
investigation  was  written  for  unstructured  grids.  The 
cells  for  all  calculations  discussed  in  this  paper  were 
generated  with  Watson’s  algorithm^  as  a  2-D  Delaunay 
triangulation  of  points  distributed  in  the  domain  and  on 
the  boundaries.  Particle  movement  and  current-cell 
identification  were  performed  simultaneously  with  a  tra¬ 
jectory-tracing  scheme,  similar  to  that  proposed  by 
Dietrich,*^  where  a  particle  is  moved  along  its  trajectory 
until  it  contacts  a  cell  boundary.  It  is  then  passed  to 
another  cell,  allowed  to  leave  the  calculation,  or 
reflected  from  a  solid  surface,  as  appropriate.  Particle 


information  was  stored  locally  in  each  cell,  following 
Dietrich  and  Boyd,-"’  to  take  advantage  of  the  cache  hard¬ 
ware  of  workstation-level  computers  and  to  prepare  for 
a  planned  transition  to  parallel  machines.  Bird’s  No 
Time  Counter  (NTC)  scheme^  was  used  for  collision 
pair  selection,  with  collision  cross-sections  calculated 
according  to  his  Variable  Hard  Sphere  (VHS)  mode!.^ 

Solid  wall  incursions  were  all  treated  as  diffuse 
reflections  with  full  thermal  and  momentum  accommo¬ 
dation.  Impinging  molecules  were  therefore  re-emitted 
with  a  velocity  randomly  selected  from  a  Maxwellian 
distribution  at  the  wall  temperature  without  regard  to 
their  incoming  state. 

Inflow/outflow  faces  required  considerably  more 
effort,  particularly  for  low-speed  cases  such  as  those 
presented  in  sections  III  A  and  IIIB.  These  faces  were 
involved  in  two  stages  of  the  calculation:  particle  move¬ 
ment  and  boundary  enforcement.  During  the  movement 
stages,  particles  were  simply  removed  from  the  calcula¬ 
tion  if  they  encountered  an  I/O  face.  Because  particle 
information  was  stored  locally  and  passed  from  cell  to 
cell  with  the  particle,  this  was  a  straightforward  opera¬ 
tion;  the  particle  was  simply  passed  to  ‘nowhere’.  Dur¬ 
ing  the  boundary  enforcement  stage,  appropriate 
particles  were  introduced  at  I/O  faces  to  maintain  the 
specified  macroscopic  flow  parameters.  Pressure  was 
specified  at  both  inflow  and  outflow  edges  and  tempera¬ 
ture  and  transverse  speed  were  specified  at  the  inflow 
edge, 

DSMC’s  statistical  nature  complicates  the  boundary 
enforcement  process,  however.  To  choose  the  quantity 
and  velocity  of  the  incoming  particles,  the  non-spec ifted 
macroscopic  quantities  must  be  determined  from  the 
current  flow  state.  Unfortunately,  there  may  be  as  few  as 
20  particles  in  a  given  cell  at  a  given  time,  so  sampling 
the  cell  containing  the  I/O  face  to  determine  this  state 
yields  a  poor  estimate.  This  leaves  two  options:  neigh¬ 
boring  cells  may  be  included  in  the  sample  or  an  aver¬ 
age  over  some  length  of  time  may  be  substituted  for  the 
instantaneous  average.  The  former  option  involves  the 
determination  of  which  neighboring  cells  have  states 
that  are  ‘close-enough’  to  that  of  the  cell  in  question  to 
yield  a  meaningful  spatial  average.  The  latter  option 
involves  choosing  an  averaging  method  that  yields  a 
sufficiently  accurate  estimate  but  still  allows  the  flow  to 
properly  reach  its  steady  state  in  a  reasonable  time 
period. 

The  latter  option  was  implemented  in  the  current 
code.  The  cell  state  was  sampled  after  a  movement  step 
was  completed,  incoming  particles  were  introduced  at 
boundary  enforcement,  and  collisions  were  performed. 
A  weighted  average  was  then  taken  between  this  result 
and  a  running  value.  The  weighting  of  the  instantaneous 
state  may  be  varied  to  balance  the  speed  of  convergence 
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with  the  accuracy  of  the  running  estimate,  A  weight  of 
1/20  was  chosen  for  the  cases  presented  in  this  paper. 

All  quantities  are  nondimensionalized  in  the  code 
and  its  output.  This  generalizes  the  construction  of  cal¬ 
culations  and  facilitates  the  interpretation  of  their 
results.  Due  to  the  particulate  nature  of  the  scheme  and 
the  importance  of  the  Knudsen  number  in  the  subject 
geometries,  the  mean  free  path  at  a  reference  condition 
is  a  natural  choice  for  nondimensionalizing  lengths.  By 
similar  reasoning,  one  of  the  molecular  speeds  (i.e, 
mean,  rms,  or  most  probable)  is  a  logical  choice  for  nor¬ 
malizing  velocities.  The  most  probable  molecular  speed 
was  selected  for  the  current  code.  Times  are  nondimen¬ 
sionalized  by  the  quotient  of  these  factors.  Number  den¬ 
sities,  pressures,  and  temperatures  are  each 
nondimensionalized  by  their  value  at  the  reference  con¬ 
dition. 

III.  Results 


reflection),  to  1,  for  full  accommodation  (diffuse 
reflection). 

Through  this  analysis,  an  expression  may  be 
obtained  for  the  pressure  distribution  in  a  microchannel 
with  diffusely-reflecting  walls  as  a  function  of  stream- 
wise  channel  location  and  overall  pressure  ratio: 

T{x)  =  (2) 

a/(  6  J  -  Z  [(  ] 

where  ^x)  and  are  the  local  and  inlet  pressures, 
respectively,  normalized  by  the  inlet  value,  Kn^,  is  the 
outlet  Knudsen  number,  x  is  the  streamwise  coordinate, 
and  L  is  the  channel  length. 

The  distribution  predicted  by  Eq.  2  may  be  com¬ 
pared  to  a  DSMC  result  as  a  test  of  both  theory  and 
code.  Such  a  comparison  is  presented  in  Fig.  i  for  a 
900x30 X  channel  with  a  pressure  ratio  of  2.4. 


A.  Slip  Flow  Regime  Micro-Channel 

The  first  case  explored  was  a  steady  flow  through  a 
micro-channel  with  an  outlet  Knudsen  number  of  0.05, 
placing  it  in  the  slip  flow  regime.  This  geometry  is  simi¬ 
lar  to  those  investigated  experimentally  by  Harley  et  al} 
and  Arkilic  et  al.^  and  numerically  (with  spectral  ele¬ 
ment  methods)  by  Beskok  and  Karniadakis.*®  This  is, 
historically,  an  important  canonical  case  for  determining 
the  effect  of  rarefaction  on  the  transport  terms  in  the 
Navier-Stokes  equations.  It  is  also  useful  as  a  represen¬ 
tation  of  the  flow  along  certain  features  common  in 
MEMS  devices,  such  as  the  space  under  the  floating 
plate  of  a  shear  stress  sensor  or  accelerometer,  which  is 
typically  only  one  micron  high  but  hundreds  of  microns 
in  breadth  and  depth.  In  addition,  it  serves  as  an  interest¬ 
ing  calibration  case  to  assess  the  accuracy  of  the  numer¬ 
ical  algorithm  because  analytical  solutions  have  been 
developed  for  this  geometry. 

In  one  such  effort,  Arkilic  et  ai  show  that  the 
Navier-Stokes  equations  may  be  solved  analytically  for 
a  long,  high  aspect-ratio,  isothermal  channel  in  the  slip 
flow  regime  if  the  boundary  conditions  are  modified  to 
include  a  Kn-dependent  streamwise  velocity  (slip)  at  the 
wall,  given  by: 


Fig.  1:  Comparison  of  computed  and  analytical  pressure 
distributions  for  a  micro-channel  in  the  slip  flow  regime. 
Excellent  agreement  is  obtained  between  code  and  theory 
Including  the  non-linear  pressure  distribution  resulting  from 
compressibility. 

A  good  agreement  is  observed,  including  the  nonlinear 
pressure  distribution  that  occurs  due  to  the  large 
pressure  drop  down  the  length  of  the  channel. 

A  theoretical  expression  for  the  streamwise  velocity 
distribution  was  also  developed  by  Arkilic  etal.: 


^walt 


2-F 

F 


wall 


(1) 


where  u  is  the  streamwise  velocity,  Kn  is  the  local 
Knudsen  number,  and  y  is  the  transverse  coordinate, 
which  has  its  zero  at  the  channel  centerline.  F  is  the 
tangential  momentum  accommodation  coefficient  and 
varies  from  0,  for  no  accommodation  (specular 


where  |i  is  the  coefficient  of  viscosity,  p  is  the  pressure, 
and  H  is  the  channel  height. 

This  equation  is  plotted  in  Fig.  2  for  the  geometry 
and  conditions  used  above.  Several  features  unique  to  a 
flow  of  this  type  are  visible.  First,  the  fluid  accelerates 
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as  it  moves  down  the  channel,  unlike  in  the  familiar  Poi- 
seuille  result.  This  is  a  consequence  of  the  density  drop 
caused  by  the  decreasing  pressure  in  the  stream  wise 
direction  (the  flow  is  effectively  isothermal).  The  mean 
streamwise  velocity  must  therefore  increase  to  maintain 
a  constant  mass  flow.  Second,  the  velocity  at  the  walls  is 
nonzero  and  increases  with  increasing  x-coordinate. 
This  is  the  aforementioned  ‘slip  flow’,  which,  by  Eq.  1, 
is  essentially  zero  for  continuum  flows  due  to  their  very 
small  Knudsen  number.  The  increase  in  slip  velocity 
down  the  channel  is  a  result  of  growth  in  both  Kn  (from 
the  decreasing  pressure)  and  velocity  gradient  (from  the 
accelerating  flow). 


Fig.  2:  Theoretical  streamwise  velocity  distribution  for  a  micro- 
channel  in  the  slip-flow  regime,  showing  flow  acceleration  due  to 
compressibility  and  non-zero  velocity  at  the  wall  due  to 
rarefaction. 


Fig.  3:  Computed  streamwise  velocity  distribution  for  the  micro- 
channel  of  Fig.  2,  showing  good  qualitative  and  quantitative 
agreement  between  analytical  and  DSMC  results 

The  DSMC  result  for  this  configuration  is  presented 
in  Fig.  3.  Comparing  this  to  the  previous  figure,  it  may 
be  concluded  that  the  DSMC  calculation  qualitatively 
reproduces  the  mean  flow  acceleration  and  the  increas¬ 


ing  slip  flow  predicted  by  the  theory.  In  addition,  good 
quantitative  agreement  is  obtained  between  the  velocity 
distributions. 

A  further  comparison  with  the  theoretical  analysis  of 
Arkilic  et  al.  may  be  found  by  normalizing  the  velocity 
distribution  of  Eq.  3  by  the  average  velocity  at  a  given 
x-location,  obtained  by  integrating  Eq.  3  from  the  lower 
wall  to  the  upper  wall: 


dP 


(I  +6Kn) 


\2\i  dx 

Rearranging  the  resulting  expression  yields: 
-Kn  + 


(4) 


(5) 


The  left  side  of  this  equation,  which  will  be  referred 
to  as  the  “similarity  speed”,  is  a  function  of  x  and  y, 
while  the  right  is  a  function  of  y  only.  Consequently,  if 
the  slip-flow  analysis  holds,  calculating  the  similarity 
speed  using  the  local  Kn{\)  and  M(x,y)  will  yield  identi¬ 
cal  parabolas  at  all  x-stations  down  the  length  of  the 
channel. 

This  assertion  was  tested  by  computing  a  similarity 
speed  distribution  from  the  DSMC  output.  The  maxi¬ 
mum  and  minimum  of  the  result  at  each  x-location  is 
shown  in  Fig.  4.  It  should  be  noted  that  the  upper  theo¬ 
retical  line  is  not  placed  at  0.25  because  an  even  number 
of  cells  was  used  in  the  DSMC  run,  so  there  is  no  data 
point  in  the  center  of  the  channel.  Also,  the  lower  line  is 
not  at  zero  because  the  macroscopic  quantities  for  a  cell 
are  assumed  to  be  associated  with  the  cell  centroid,  so 
there  are  no  data  points  on  the  walls  themselves. 


0  5  10  15  20  25  30 
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Fig.  4:  Comparison  of  computed  and  theoretical  maximum  and 
minimum  similarity  speeds  for  a  micro-channel  in  the  slip  flow 
regime,  showing  that  similarity  is  indeed  achieved  and  that 
observed  magnitudes  compare  well  with  predicted  values. 

It  may  be  concluded  from  this  figure  that  the  similar¬ 
ity  assertion  indeed  holds  for  the  slip  flow  channel.  As 
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predicted  by  the  analysis,  the  down-channel  variation  of 
the  streamwise  velocity  profile  seen  in  Fig.  3  has  given 
way  to  a  constant  similarity  speed  profile.  In  addition, 
the  maximum  and  minimum  similarity  speeds  compare 
well  with  the  predicted  values. 

Overall,  excellent  agreement  was  obtained  between 
the  analytical  solution  of  Arkilic  et  at.  and  the  DSMC 
results.  This  supports  the  accuracy  of  both  techniques. 
For  the  DSMC  code,  however,  it  is  just  the  beginning. 
Many  more  interesting  flows,  for  which  there  are  no 
reliable  analytical  solutions,  may  be  easily  treated  with 
this  method.  The  remaining  cases  presented  in  this  paper 
are  intended  to  demonstrate  this  capability. 

B.  Transition  Regime  Micro-Channel 

One  of  the  great  strengths  of  DSMC  is  its  validity 
for  dilute  gases  in  all  Knudsen  number  regimes.  One  of 
the  most  interesting  of  these  is  the  transition  regime, 
defined  previously  as  0. 1  <  Kn  <  3.  Here  the  mean  free 
path  is  comparable  to  the  characteristic  dimension  of  the 
flow.  This  makes  analytical  solution  very  difficult 
because  the  approximation  of  transport  terms  based  on 
macroscopic  quantities  becomes  inaccurate,  precluding 
the  use  of  the  Navier-Stokes  equations  (even  with  the 
slip-flow  boundary  condition  of  Eq.  1).  Collisions  are 
still  important,  however,  so  the  collisionless  Boltzmann 
equation  is  not  yet  an  option.  This  leaves  only  the  full 
Boltzmann  equation;  a  very  difficult  expression  to  solve, 
either  analytically  or  with  standard  numerical  tech¬ 
niques. 

DSMC  is  therefore  a  very  attractive  tool  for  investi¬ 
gating  the  transition  regime.  It  is  also  one  of  the  few 
tractable  techniques  which  is  uniformly  valid  in  mixed 
Kn-regimes.  These  are  important  features  because,  due 
to  the  aforementioned  difficulties,  relatively  little  is 
known  about  these  flows.  Such  knowledge  is  critical, 
however,  because  many  MEMS  devices  contain  regions 
in  both  the  slip-flow  and  transition  regimes. 

The  micro-channel  was  again  used  as  sample  case, 
this  time  to  observe  the  failure  of  the  slip-flow  analysis 
as  the  flow  enters  the  transition  regime.  The  channel 
treated  here  is  1 20x2  A,,  with  a  pressure  ratio  of  4.2  and 
an  outlet  Knudsen  number  of  0.44.  The  working  fluid  is, 
again.  Nitrogen. 

Proceeding  as  for  the  slip-flow  case,  the  computed 
pressure  distribution  is  compared  to  the  slip-flow  pre¬ 
diction  in  Fig.  5.  The  continuum  curve  (Kn  =  0.0)  is  also 
shown  for  reference.  As  expected,  the  excellent  agree¬ 
ment  between  DSMC  and  theory  obtained  for  the  slip- 
flow  case  (Fig.  1)  is  no  longer  present.  The  error 
between  the  curves  has  grown  from  less  than  2%  to 
more  than  4%.  The  form  of  this  disagreement  is  also  sig¬ 
nificant:  the  computed  curve  is  more  linear  than  its  ana¬ 


lytical  counterpart.  A  trend  of  increasing  pressure  curve 
linearity  with  increasing  rarefaction  is  therefore  estab¬ 
lished  by  the  relative  shape  of  the  continuum,  slip-flow, 
and  transition  curves. 


Fig.  5:  Comparison  of  computed,  slip-flow  and  continuum 
pressure  distributions  for  a  micro-channel  in  the  transition 
regime.  The  transition- regime  curve  is  more  linear  than  the  slip 
flow  curve  which  is,  in  turn,  more  linear  than  the  continuum 
curve. 

The  analytical  prediction  for  the  streamwise  velocity 
distribution  in  this  channel  is  presented  in  Fig.  6.  Note 
that  the  theory  predicts  a  flatter  profile  than  for  the  pre¬ 
vious  case.  This  may  be  attributed  to  the  last  term  in  Eq. 
3,  which  is  constant  across  the  channel  and  proportional 
to  Kn.  The  Knudsen  number  is  an  order  of  magnitude 
larger  in  this  case,  so  this  term  has  a  much  stronger 
influence  upon  the  distribution. 


Fig.  6:  Theoretical  streamwise  velocity  distribution  for  a 
transition-regime  micro  channel  calculated  with  the  slip-flow 
analysis. 

Upon  plotting  the  computed  distribution  for  compar¬ 
ison  in  Fig.  7.,  it  becomes  evident  that  the  assumptions 
supporting  Eq.  3  are  beginning  to  fail.  Both  the  slip  flow 
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and  maximum  speeds  at  a  given  x-location  are  higher 
than  predicted  by  as  much  as  40%.  This  allows  the  chan¬ 
nel  to  support  a  mass  flow  greater  than  would  be  pre¬ 
dicted  by  the  slip  flow  theory,  which,  in  turn,  predicts  a 
mass  flow  greater  than  the  traditional  Navier-Stokes 
analysis. 


Fig.  7:  Computed  streamwise  velocity  distribution  for  the  channel 
of  Fig.  6.  The  magnitude  of  both  the  slip  and  the  local  maximum 
velocity  are  greater  than  the  analytical  prediction  by  as  much  as 
40%. 

A  final  exposition  of  transition  behavior  may  be 
made  via  the  similarity  analysis  of  the  previous  section. 
As  mentioned  previously,  the  transition  regime  is  com¬ 
monly  considered  to  begin  at  Kn  =  0, 1.  This  assertion 
may  be  tested  by  noting  that  Kn  increases  with  down¬ 
stream  position.  The  similarity  profiles,  which  depend 
on  the  slip-flow  solution,  may  therefore  be  computed  for 
each  position  and  compared  to  the  analytical  slip-flow 
prediction.  Because  each  position  has  a  Knudsen  num¬ 
ber  associated  with  it,  the  value  for  which  the  slip  flow 
analysis  fails  may  be  determined  by  finding  the  point 
where  the  experimental  and  analytical  curves  begin  to 
diverge.  Toward  this  end,  Fig.  8  contains  the  computed 
maximum  and  minimum  similarity  speeds,  plotted  with 
the  analytical  prediction  in  a  fashion  identical  to  that  of 
Fig.  4.  The  Kn  distribution  was  then  overlaid  to  facili¬ 
tate  determining  its  value  when  the  slip-flow  analysis 
fails. 

It  may  be  seen  in  this  figure  that  the  slip  flow  analy¬ 
sis  begins  to  fail  at  approximately  Kn  =  0.15.  This  sup¬ 
ports  the  oft-used  boundary  for  the  transition  region,  Kn 
=  0.1.  This  limit  may  be  understood  if  the  wall  boundary 
condition,  Eq.  1,  is  viewed  as  an  expansion  of  the  wall 
velocity  in  powers  of  Kn.  The  no-slip  condition  is  then 
the  zeroth-order  solution,  and  the  slip  condition  of  Eq.  I 
is  the  first-order  solution.  It  is  therefore  logical  that  the 
neglected  higher-order  terms  would  begin  to  signifi¬ 
cantly  affect  the  result  when  Kn  exceeds  0. 1.  A  second- 
order  accurate  boundary  condition  in  terms  of  the  con¬ 


tinuum  variables  is  presented  in  Ref  10,  though  it  should 
be  remembered  that  the  Navier-Stokes  equations  them¬ 
selves  are  only  strictly  valid  to  first  order  in  Kn. 


-  DSMC 

-  -  -  -  Prediction  Kn  — — — 


Fig.  8:  Computed  and  analytical  similarity  speeds  with  Knudsen 
number  overlay.  The  point  at  which  the  computed  curve  begins 
to  diverge  from  the  prediction  corresponds  well  with  the 
commonly-stated  Kn  for  the  onset  of  the  transition  regime  (0.1). 

C.  Supersonic  Micro-Nozzle 

The  final  case  presented  is  a  supersonic  micro- noz¬ 
zle.  This  is  a  channel  whose  upper  and  lower  walls  form 
a  parabolic  contraction/expansion.  The  grid  for  this  cal¬ 
culation  is  shown  in  Fig.  9. 


0  1  2  3  4  5  6 

x/HJ 

Fig.  9:  Computational  grid  for  supersonic  micro-nozzle.  The 
upper  and  lower  walls  are  parabolic  and  the  area  ratio  is  3.5. 

For  a  sufficiently  long  nozzle,  this  geometry  would 
be  quasi- ID  and  the  analysis  used  in  section  IIIA,  with 
appropriate  modifications  for  the  slowly  varying  chan¬ 
nel  height,  would  be  valid.  For  the  case  shown  in  Fig.  9, 
however,  the  total  length  of  the  nozzle  is  only  six  times 
the  throat  height,  so  the  quasi- ID  assumption  is  not 
applicable.  In  addition,  the  significant  expansion  of  the 
fluid  may  also  cause  the  Kn-regime  to  change  at  one  or 
more  streamwise  positions.  These  factors  make  analyti¬ 
cal  and  continuum-based  numerical  treatment  of  this 
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geometry  difficult.  Nozzles  such  as  these  may  play 
important  roles  in  micro  devices  like  micro-rocket 
thrusters  and  micro-gas  turbine  generators,  for  example. 
Investigating  their  behavior  is  therefore  a  valuable  task 
for  which  DSMC  is  well-suited. 

This  particular  nozzle  was  run  with  an  inlet  pressure 
of  10  times  the  reference  pressure  and  exhausted  to 
‘vacuum’.  The  latter  condition  was  implemented  by 
simply  removing  from  the  simulation  any  particle  which 
crossed  the  outlet  boundary  and  refraining  from  intro¬ 
ducing  new  particles  on  those  faces.  The  resulting  pres¬ 
sure  ratio  was  approximately  24  and  the  outlet  Knudsen 
number,  based  on  the  passage  height,  was  0.03.  The 
walls  were  isothermal  at  the  reference  temperature  with 
full  energy  and  momentum  accommodation.  The  work¬ 
ing  fluid  was  Helium,  which  was  supplied  at  the  refer¬ 
ence  temperature. 

The  Mach  number  distribution  for  this  case  is  shown 
in  Fig.  10.  A  number  of  interesting  features  are  visible. 
First,  as  expected,  the  flow  is  sonic  at  the  throat.  Second, 
a  Mach  number  of  2.4  is  reached.  This  is  considerable 
when  it  is  noted  that  the  nozzle  is  only  approximately 
600  inlet  mean  free  paths  in  length.  Finally,  the  slip  flow 
speed  is  substantial,  exceeding  Mach  0.5  near  the  outlet. 


Fig.  10:  Mach  number  distribution  for  a  supersonic  micro- nozzle. 
A  maximum  Mach  number  of  2.4  is  obtained  from  an  inlet 
pressure  of  10  and  a  vacuum  outlet  The  slip  velocity  exceeds 
Mach  0.5  near  the  outlet 

The  temperature  distribution  for  the  micro  nozzle  is 
shown  in  Fig.  1 1 .  It  is  clear  that,  unlike  the  channel  case, 
this  flow  cannot  be  considered  isothermal.  A  strong  tem¬ 
perature  gradient  exists  in  both  the  stream  wise  and 
transverse  directions,  creating  another  obstacle  to  ana¬ 
lytical  treatment.  This  is  a  very  notable  feature  when  the 
diminutive  dimensions  of  the  nozzle  are  considered. 
Though  the  walls  are  only  ten  mean  free  paths  apart  at 
the  throat  and  are  isothermal  with  full  energy  accommo¬ 
dation,  the  fluid  is  still  able  to  realize  a  substantial 


reduction  in  temperature.  The  effect  of  rarefaction  has 
therefore  been  to  significantly  reduce  the  thermal  com¬ 
munication  between  the  wall  and  the  fluid.  This  asser¬ 
tion  is  supported  by  noting  the  large  thermal  slip  at  the 
wall. 


Fig.  11:  Temperature  distribution  for  a  supersonic  micro-nozzle.  A 
significant  temperature  jump  is  present  despite  the  fact  that  the 
walls  are  isothermal  (at  temperature  1.0)  with  full  thermal 
accommodation. 

IV.  Conclusion 

The  direct  simulation  Monte  Carlo  method  has 
proven  to  be  a  valuable  tool  for  investigating  the  behav^ 
ior  of  flows  which  are  considered  ‘rarefied’  due  to  min¬ 
iaturization.  This  paper  has  presented  only  a  small 
subset  of  the  many  cases  of  interest  to  both  the  fluid- 
dynamicist  and  the  MEMS  designer  which  are  easily 
treated  with  this  method. 

In  section  IIIA,  DSMC  results  for  a  slip-flow  regime 
micro-channel  were  compared  to  an  analytical  solution 
of  the  Navier-Stokes  equations  presented  by  Arkilic  et 
al.  Several  aspects  of  computed  and  theoretical  stream- 
wise  velocity  and  pressure  distributions  were  consid¬ 
ered.  In  all  cases,  the  numerical  results  were  found  to 
compare  well  with  their  analytical  predictions. 

By  manipulating  the  analytical  results,  an  expression 
for  a  streamwise  velocity  ‘similarity’  profile  was  also 
developed.  It  was  subsequently  found  that  the  computed 
velocity  distribution,  which  is  a  function  of  x,  collapsed 
quite  well  to  this  x-independent  profile  when  the  flow 
variables  were  processed  according  to  the  theoretical 
expression. 

These  findings  support  the  accuracy  of  both  the  code 
and  analytical  work  and  demonstrate  another  valuable 
function  of  DSMC  computations:  the  evaluation  of  ana¬ 
lytical  models  for  rarefied  flow  behavior.  Developing 
these  models  is  an  important  task  because  a  great  num¬ 
ber  of  MEMS  geometries  fall  into  the  slip-flow  regime. 
Analytical  tools  have  great  value  to  the  designer 
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because  they  are  much  less  expensive  than  a  full  DSMC 
calculation  but  are  able  to  approximate  behavior  invisi¬ 
ble  to  continuum  techniques.  They  are  difficult  to  vali¬ 
date,  however,  because  the  effects  they  describe  are 
often  too  small  for  effective  experimental  investigation. 

In  section  IIIB,  a  DSMC  computation  for  a  transition 
regime  micro-channel  was  compared  to  predictions 
from  the  slip-flow  expressions  of  Arkilic  et  ai  By  con¬ 
sidering  the  form  of  disagreement  between  these  results, 
the  behavior  of  transition-regime  flows  was  illuminated. 
This  is  an  especially  important  application  of  DSMC 
because  very  few  techniques,  either  analytical  or  numer¬ 
ical,  are  available  for  analyzing  these  flows.  This  is  a 
significant  problem  for  the  MEMS  community  because 
the  geometries  and  working  conditions  for  many 
devices  place  them  in  this  regime. 

In  section  IIIC,  a  parabolic  micro-nozzle  was  inves¬ 
tigated.  This  case  is  intended  to  represent  devices  which 
contain  complexities  not  amenable  to  other  forms  of 
solution.  Its  more  complicated  geometry,  lack  of  isother¬ 
mal  flow,  and  sharp  gradients  pose  problems  for  many 
types  of  analysis.  No  special  considerations  for  these 
features  were  required  to  perform  the  DSMC  calcula¬ 
tion,  however.  It  is  therefore  demonstrated  to  be  a  very 
versatile  method  whose  value  increases  with  flow  com¬ 
plexity. 

The  unique  features  of  rarefied  gas  flows  have  many 
implications  for  MEMS  designers.  Their  larger  mass 
flow  rate  for  a  given  geometry  and  inlet  conditions  must 
be  considered  when  designing  control  systems  for 
micro-chemical  reactors  and  the  decreased  thermal  com¬ 
munication  between  the  fluid  and  its  surroundings  must 
be  considered  for  applications  involving  temperature 
measuring  devices  or  heaters,  for  example. 

DSMC  therefore  has  great  promise  for  improving 
the  design  and  optimization  of  MEMS.  Its  unique 
strength  for  this  field  is  its  ability  to  calculate  in  any  of 
the  four  Knudsen  number  regimes  without  modification. 
This  is  especially  valuable  for  flows  that  have  regions  in 
different  regimes,  as  is  the  case  for  some  MEMS,  The 
addition  of  unstructured  grid  capability  and  the  trajec¬ 
tory-tracing  scheme  give  the  code  the  ability  to  handle 
arbitrary  geometries.  This  is  a  valuable  asset  when  ana¬ 
lyzing  the  complex  structures  included  in  many  of  these 
devices.  In  addition,  it  is  quite  straightforward  to 
include  flow  features  such  as  chemical  reactions,  multi¬ 
species  mixing  and  particle  transport  in  the  code  due  to 
its  modular  nature.  Finally,  its  cell-based  structure 
makes  it  well-suited  for  parallel  computation,  which  is 
an  increasingly-important  attribute  for  a  large-scale 
numerical  scheme. 
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Abstract 

Numerical  results  on  the  convergence  of  the  Di¬ 
rect  Simulation  Monte  Carlo  (DSMC)  technique  are 
presented  for  both  well-resolved  cases  and  cases  in 
which  the  cell  size  is  much  larger  than  the  mean- 
free-path  of  the  gas.  Results  from  two  test  cases 
are  presented:  a  one-dimensional  Rayleigh  flow  and 
a  two-dimensional  spatially-developing  shear  layer. 
Two  main  effects  are  observed  as  the  cell  Knudsen 
number  decreases.  Firstly,  the  particle  population 
in  each  cell  approaches  an  equilibrium  distribution, 
suggesting  that  the  DSMC  is  formally  solving  the 
Euler  equations,  not  the  Navier-Stokes  equations. 
Secondly,  an  artificial  viscosity  whose  magnitude  is 
inversely  proportional  to  the  cell  Knudsen  number  is 
observed  to  produce  unreasonably  thick  shear  lay¬ 
ers.  This  viscosity  is  due  to  thermal  fluctuations 
which  cause  an  unphysical  diffusion  of  momentum 
across  cell  boundaries  as  the  cell  size  increases.  This 
diffusion  (related  to  the  position-independent  na¬ 
ture  of  the  DSMC  collision  algorithm)  decreases  at 
higher  Mach  numbers  as  the  relative  magnitude  of 
the  thermal  fluctuations  decreases.  Numerical  re¬ 
sults  are  presented  showing  the  ability  of  the  DSMC 
method  to  accurately  solve  complex  (inviscid)  con¬ 
tinuum  flows. 

Nomenclature 

c  Most-probable  molecular  speed 

/  Velocity  distribution  function 

H  .  Computational  domain  height 
fo  Reference  (Maxwellian)  distribution 

K  Rayleigh-problem  Knudsen  number: 

X/y/u7rt 

Kc  Knudsen  Number  based  on  cell:  A/Az 

Ma  Mach  number 
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Np  Number  of  Collision  pairs 

Nc  Number  of  particles-per-cell 

Nref  Reference  number  density 

R  Gas  Constant 

Re  Reynolds  number 

t  Time 

T  Temperature 

U  Streamwise  velocity 

Streamwise  and  wall-normal  perturbation 
velocities 

Uo  Wall  Speed 

Ui,U2  Upper  and  Lower  speeds  in  shear  layer 

Vc  Cell  Volume  (normalized  by  A) 

x,2/  Streamwise  and  wall-normal  coordinates 

Axd  Dimensional  cell  size 

Ax  Non-dimensional  cell  size:  Ax^/A  =  1/Kc 
Atd  Dimensional  time-step 

A^  Non-dimensional  time  step:  cAi/ A 

Tj  Similarity  coordinate:  y/{2y/iH) 

A  Mean  free  path  of  gas 

jy  Kinematic  viscosity 

6  Momentum  thickness 

1  Introduction 

Thanks  to  numerous  advances  in  Computational 
Fluid  Dynamics  (CFD)  over  the  past  thirty  years, 
it  is  now  relatively  straightforward  to  compute 
problems  of  engineering  interest  in  two  and  three- 
dimensions,  solving  both  the  Euler  equations  and, 
to  a  lesser  extent,  the  Navier-Stokes  equations.  In 
parallel  to  this  development,  numerical  techniques 
for  rarefied  gas  flows  have  also  developed  rapidly  in 
the  past  several  years.  Although  a  variety  of  tech¬ 
niques  have  been  used,  including  Molecular  Dynam¬ 
ics  [7]  and  solution  techniques  based  on  the  Bur¬ 
nett  equations  [8],  the  most  popular  has  been  the 
Direct  Simulation  Monte  Carlo  (DSMC)  technique, 
popularized  primarily  by  Bird  [3].  In  the  DSMC 
approach,  the  evolution  of  test  particles,  each  rep¬ 
resentative  of  a  large  number  of  real  molecules,  is 
computed.  While  test  particle  motion  is  computed 
exactly,  collisions  between  test  particles  are  treated 


in  a  probabilistic  manner.  Although  the  technique 
is  not  formally  identical  to  the  exact  gas  behavior,  it 
has  been  shown  to  be  useful  for  problems  spanning 
a  wide  range  of  flow  conditions  ranging  from  one- 
,  dimensional  shock  structure  [3]  to  the  flow  over  the 
Magellan  Spacecraft  as  it  descends  into  the  Venusian 
atmosphere  [5].  Although  the  method  is  motivated 
by  problems  in  rarefied  flow  regimes,  it  is  of  interest 
to  examine  how  the  method  handles  flows  with  vam 
ishing  Knudsen  number,  K,  that  is,  flows  in  which 
the  characteristic  geometry  is  much  greater  than  the 
mean  free  path  of  the  working  gas.  The  reasons  for 
this  interest  is  for  computations  which  contain  both 
continuum  and  rarefied  regimes  as  well  as  for  use 
in  fully  continuum  flows  where  a  particle  approach 
might  offer  some  advantages  over  a  field  equation 
formulation.  It  is  this  issue  that  is  discussed  in  the 
current  paper. 


2  Scaling  issues  for  DSMC 

Central  to  the  DSMC  approximation  is  the  as¬ 
sumption  that  gas  properties,  such  as  temperature, 
pressure,  etc.  are  constant  across  each  cell.  This 
allows  for  collision  pairs  selected  within  a  cell  to  be 
independent  of  position  in  the  cell  and  requires  that 
the  cell  size  be  smaller  than  the  gradients  in  the 
flow  field.  This  is  the  first  condition  that  must  be 
satisfied  in  scaling  to  low-Knudsen  number  flows  and 
can  be  satisfied  without  any  difficulty.  An  example 
of  what  happens  when  this  is  violated  is  discussed 
in  the  section  4.1  with  regard  to  the  resolution  of 
shock  thicknesses.  A  more  subtle  issue,  however,  is 
the  scaling  of  the  time  step  and  particle  collisions 
and  this  is  discussed  in  the  following  section. 


this  constraint  may  be  conveniently  expressed  as  a 
“CFL”,  or  Courant-Friedrich-Lewy  condition: 


cAtd 

Axd 


<  1. 


(3) 


The  basis  for  this  constraint  is  that  a  particle  should 
reside  in  a  given  cell  for  a  few  time-steps  so  that 
there  are  enough  opportunities  for  each  particle  to 
effectively  interact  with  other  particles  so  that  in¬ 
formation  can  be  distributed  throughout  the  com¬ 
putational  domain.  It  should  be  realized  that,  un¬ 
like  its  CFD  counterpart,  the  CFL  condition  for 
DSMC  computations  is  not  a  stability  requirement, 
but  rather  a  validity  requirement.  Violation  of  the 
condition  will  still  yield  results,  but  they  may  be 
inaccurate  due  to  the  lack  of  residence  time  for  a 
particle  in  each  cell. 


Inherent  in  the  DSMC  method  is  the  decoupling  of 
particle  movement  and  their  interactions  (collisions, 
reactions,  etc.).  Once  the  particles  have  moved  their 
collisions  may  be  computed  though  a  variety  of  tech¬ 
niques.  A  common  technique  is  the  so-called  No- 
Time-Counter  (NTC)  method  [2],  in  which  the  num¬ 
ber  of  collision  pairs,  iVp,  to  be  considered  in  a  given 
cell  is  computed  through  the  equation: 


Np  oc  At 


T^ref^c 


(4) 


where  K  is  the  volume  of  a  cell,  Nc  is  the  number 
of  test  particles  in  the  cell,  Uref  is  the  simulation 
number  density  (the  total  number  of  simulated  par¬ 
ticles  divided  by  the  total  simulation  volume).  All 
quantities  are  now  normalized  by  A  and  c.  By  realiz¬ 
ing  that  the  normalized  cell  volume  scales  like  K~^ 
and  that  the  simulation  number  density  scales  like 
NcK^,,  we  see  that 


A  typical  rule-of-thumb  for  the  DSMC  technique 
is  that  the  typical  cell  size  be  of  the  order  of  the 
mean-free-pat h  of  the  gas,  A.  This  can  be  expressed 
as  a  condition  on  the  “cell  Knudsen  number”,  Kc'. 


Kc  = 


(1) 


where  Axd  is  a  typical  cell  size  dimension  (the  sub¬ 
script  “d”  referring  to  dimensional  quantities.  A 
second  requirement  is  that  the  time  step  be  small 
compared  to  a  characteristic  time: 


Atd  < 


Axd 


(2) 


^  «  At  (5) 

In  other  words,  for  a  fixed  size  computation  [Nc  held 
constant),  Np  depends  only  on  the  non-dimensional 
time  step  chosen  to  advance  the  solution.  We  term 
this  ratio  the  “over-collision  ratio”  since  it  is  an  indi¬ 
cation  of  the  ratio  of  collision  pairs  to  the  number  of 
particles  in  the  cell  and  a  value  equal  to  one  implies 
that  every  particle  will  considered  for  a  collision  in¬ 
teraction  during  each  time  step.  In  a  well-resolved 
DSMC  computation,  where  At  «  0.2,  approximately 
20%  of  a  given’s  cell’s  particles  will  be  considered  for 
collisions. 


where  c  is  the  most  probable  molecular  speed.  Bor¬ 
rowing  from  traditional  computational  mechanics, 
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The  over-collision  ratio  reveals  the  essential  scal¬ 
ing  problem  when  applying  DSMC  methods  to  low- 
Knudsen  number  flows.  As  the  size  of  the  problem 


increases,  it  becomes  computationally  impractical  to 
continue  to  maintain  a  cell  size  that  is  of  the  order 
of  the  mean-free-path.  Given  that,  one  is  forced  to 
let  the  cell  size  increase  or,  in  the  framework  out¬ 
lined  above,  to  let  Kc  decrease  below  one.  As  the 
cell  size  increases,  the  size  of  the  time  step  needs 
to  be  reconsidered.  We  have  two  options:  By  scal¬ 
ing  the  time  step  so  that  the  CFL  number  remains 
constant,  information  propagation  across  the  com¬ 
putational  domain  and  computational  efficiency  is 
preserved.  However,  this  option  implies  that  the 
number  of  collision  pairs  to  be  considered  quickly  be¬ 
comes  extremely  large,  resulting  in  an  over-collision 
ratio  larger  than  one  (in  fact,  much  larger)  as  well  as 
a  computationally  unreasonable  number  of  collisions 
partners  to  process,  A  second  option  is  to  maintain 
Ai  at  some  small  value.  While  this  maintains  a  rea¬ 
sonable  value  of  Np/Nc^  it  results  in  a  very  small 
value  of  the  CFL  number  which  results  in  an  ineffi¬ 
cient  advancement  of  the  solution  and  accumulation 
of  statistics.  Additionally,  one  can  argue  effectively 
that  a  small  CFL  number  makes  no  sense  for  it  im¬ 
plies  that  most  particles  will  take  many  time  steps 
to  cross  a  given  cell.  Given  this,  the  collision  phase 
of  each  time  step  will  involve  the  same  group  of  par¬ 
ticles  as  in  the  previous  time  step  since  almost  no 
particles  leave  or  enter  a  cell.  The  net  effect,  there¬ 
fore,  is  as  if  a  larger  time  step  were  used,  but  with  the 
computationally  pointless  exercise  of  moving  all  the 
particles  some  small  distance  in-between  each  colli¬ 
sion  phase  of  the  calculation.  For  this  reason,  we 
argue  that  it  only  makes  sense  to  maintain  a  con¬ 
stant  CFL  number,  regardless  of  the  cell  Knudsen 
number,  Kc- 


A  solution  to  this  problem  was  proposed  by  Bartel 
et  al  [1]  who  recognized  that  the  large  number  of 
collisions  in  a  cell  suggested  by  (4)  do  not  serve  any 
purpose  other  than  to  reinforce  a  Maxwellian  distri¬ 
bution  amongst  the  particles  in  a  given  cell.  Given 
this,  it  should  be  sufficient  to  restrict  the  number 
of  collisions  to  some  small  number  (comparable  to 
the  number  of  particles  in  the  cell)  and  still  take 
large  time  steps  during  the  computation.  In  the  cur¬ 
rent  terms,  Bartel’s  suggestion  was  to  limit  the  over¬ 
collision  ratio  to  some  (arbitrary)  value  less  than  its 
“true”  value  given  by  (4).  This  approach  was  illus¬ 
trated  by  Bartel  et  al.  for  a  Couette  flow  and  an 
expanding  nozzle  flow  with  good  results.  One  pur¬ 
pose  of  this  paper  is  to  examine  the  implications  of 
the  collision  limiter  under  somewhat  closer  scrutiny. 


Uo 


Figure  1:  Schematic  of  the  One-Dimensional 
Rayleigh  flow 

3  Numerical  Results 

To  explore  these  ideas,  one-dimensional  unsteady 
Rayleigh  flow  was  chosen  as  a  model  problem.  This 
flow  (illustrated  in  Figure  1)  consists  of  a  stationary 
gas  subjected  to  the  sudden  motion,  C/o,  of  the  wall. 
The  problem  is  well-suited  for  this  study  since  it  does 
not  have  an  length  scale  imposed  by  geometry,  but 
rather  by  time.  In  addition,  it  is  a  one-dimensional 
flow  so  computations  are  quick,  allowing  for  several 
test  cases  to  be  considered. 


3.1  Analytic  Considerations 


We  can  solve  this  problem  analytically  in  two  dis¬ 
tinct  regimes.  If  the  Knudsen  number  is  high,  the 
collisionless  Boltzmann  solution  is  given  by: 


y 


V2RTt, 


(6) 


where  R  is  the  gas  constant  and  T,  the  temperature 
and  erfc  is  the  complimentary  error  function.  At  the 
other  extreme,  if  the  Knudsen  number  is  small,  the 
Navier-Stokes  equation  (for  incompressible  flow)  for 
this  geometry  is  given  by: 


du  _ 


(7) 


where  u  is  the  kinematic  viscosity  (proportional  to 
the  mean-free-path,  A).  This  equation  can  be  solved 
with  a  slip-flow  boundary  condition: 


resulting  in: 


du 

^  y=0 

(8) 

rr  erfc(77) 

^°K  +  1 

(9) 
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where  /y  is  a  similarity  variable: 


T)  = 


(10) 


and  K  is  the  Knudsen  number,  here  defined  in  a 
somewhat  unusual,  but  convenient  manner  as: 


K  = 


(11) 


This  is  a  more  general  version  of  the  classical 
Rayleigh  solution  [6].  Note  that  in  this  solution,  K 
is  a  function  of  time  and  becomes  infinite  as  t  — >•  0. 
This  makes  intuitive  sense  since  the  characteristic 
scale  of  the  solution  is  the  momentum  thickness  of 
the  viscous  layer  which  is  initially  zero,  but  grows 
as  time  increases.  One  implication  of  this  time- 
dependent  Knudsen  number,  however,  is  that  some 
care  must  be  taken  in  evaluating  the  solution  at 
77  =  0  and  00  where  appropriate  limits  of  both  t 
and  y  must  be  taken. 


3.2  Well-Resolved  Computations 

This  section  presents  DSMC  results  for  cases  that 
were  well-resolved  -  i.e.  where  the  cell  Knudsen  num¬ 
ber  was  greater  than  or  equal  to  one.  For  these 
case  the  CFL  number  used  was  0.3.  In  order  to  en¬ 
sure  that  the  incompressible  Navier-Stokes  solution 
was  valid,  the  wall  velocity  was  chosen  to  be  0.2 
{Ma  =  0.22).  The  DSMC  computations  were  car¬ 
ried  out  using  a  two-dimensional  DSMC  code  writ¬ 
ten  for  unstructured  grids.  The  gas  for  all  cases 
presented  here  is  Helium  and  the  collision  cross- 
sections  were  computed  using  the  Variable- Hard- 
Sphere  model  of  Bird  [2]. 

Figure  2  shows  the  near- wall  velocity  profile,  U (y), 
for  very  small  times  {t  =  0.004)  i.e.  times  for  which 
the  collisionless  Boltzmann  solution  should  still  be 
valid.  The  solid  line  represents  the  analytic  solution 
given  by  (6)  while  the  symbols  represent  the  DSMC 
solution.  For  this  computation,  the  cell  Knudsen 
number  was  1.3  x  10^  and  the  time-step  2.5  x  10“^. 
We  see  that  the  agreement  is  excellent. 

Figure  3  shows  the  velocity  profile  at  a  later  time 
{t  =  100)  where  we  would  expect  the  continuum  so¬ 
lution  to  apply.  Note  that  at  this  time,  K  =  0.07  and 
there  is  still  a  perceptible  velocity  slip  at  y  =  0.  The 
DSMC  results  are  now  compared  with  the  Navier- 
Stokes  solution  (9)  and  show  good  agreement. 

If,  as  figure  3  suggests,  the  DSMC  is  correctly  solv¬ 
ing  the  Navier-Stokes  equations,  the  velocity  distri¬ 
bution  should  reflect  this  as  a  slight  perturbation 


Figure  2:  Velocity  profile,  U{y)  a,t  t  =  0.004(jK  = 
11.2).  The  symbols  correspond  to  the  numerical 
solution  from  a  well-resolved  DSMC  computation 
while  the  solid  line  is  the  analytic  solution  to  the 
collisionless  Boltzmann  equation  (6).  the  y  axis  is 
normalized  by  A. 


Figure  3:  Velocity  profile,  U{y)  at  i  =  100  {K  = 
0.07).  The  symbols  correspond  to  the  numerical 
solution  from  a  well-resolved  DSMC  computation 
while  the  solid  line  is  the  analytic  solution  to  the 
slip-flow  Navier-Stokes  equations  (9). 
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Figure  4:  Distribution  of  velocity  perturbations  for 
a  well-resolved  computation,  y  =  6.14,  t  =  40,  The 
reference  (Maxwellian)  distribution  has  been  sub¬ 
tracted  and  the  four  quadrants  of  have  been 

collapsed  onto  one  using  the  known  anti-symmetry 
of  the  Chapman-Enskog  distribution. 


from  Maxwellian,  given  by  the  Chapman-Enskog 
distribution  for  a  one-dimensional  isothermal  shear 
flow  [3]: 
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Figure  5:  DSMC-computed  viscosity  as  a  function  of 
time  for  a  well-resolved  computation.  Also  plotted 
is  the  computed  Knudsen  number,  K. 


viscosity,  but  that  it  appears  to  asymptote  to  the 
correct  value  as  K  ^0.  The  variation  in  the  data, 
and  the  slight  overprediction  of  the  viscosity  at  later 
times  is  most  likely  due  to  a  variety  of  factors  includ¬ 
ing  statistical  scatter  due  to  the  low-Mach  number 
and,  at  later  times,  a  weak  influence  of  the  far  wall 
(located  200 A  from  the  moving  surface). 


f{u',v')  =  fo 


uu'v'  dU\ 
{RTr  dy) 


(12) 


where  fo  is  the  equilibrium  distribution.  This  was 
computed  and  is  shown  in  figure  4,  sampled  at 
y  =  6.14,  i  =  40.  Here,  the  average  mean  veloc¬ 
ity,  [/,  has  been  subtracted  and  the  four  quadrants 
of  {u\v')  have  been  collapsed  onto  one,  preserving 
the  appropriate  anti-symmetry  by  subtracting  across 
the  axes.  The  Chapman-Enskog  distribution  is  for¬ 
mally  only  valid  for  small  perturbations  compared  to 
unity  and  since  this  is  a  low  Mach  number  compu¬ 
tation,  the  thermal  fluctuations  are  not  that  small 
and  so  only  qualitative  comparisons  are  appropri¬ 
ate.  Nevertheless,  we  see  that  the  deviation  from 
Maxwellian  computed  in  the  DSMC  simulation  is  in 
good  agreement  with  its  predicted  structure  indicat¬ 
ing,  as  expected,  that  the  gas  is  weakly  perturbed 
from  equilibrium  by  the  shear. 


Figure  5  shows  the  time-dependence  of  the  vis¬ 
cosity  as  the  computation  evolves  from  t  =  0 
to  200.  The  viscosity  is  computed  via  a  nonlin¬ 
ear  least-squares  fit  of  the  analytic  solution  (9)  to 
the  numerically-derived  velocity  profile.  Given  the 
non-dimensionalization  described  above,  the  correct 
value  for  the  viscosity  of  Helium  is  0.64  and  we  see 
that  for  small  times,  the  solution  under-predicts  the 


The  good  agreement  between  DSMC  and  analytic 
results  presented  in  this  section  are  not  surprising 
since  the  success  of  the  DSMC  technique  in  modeling 
the  true  gas  physics  is  well-documented  (e.g.  Bird, 
[3]).  We  present  these  results  purely  for  calibration 
purposes.  What  is  of  real  interest  is  to  investigate 
what  happens  cts  the  cell  Knudsen  number  decreases 
below  one  or,  in  other  words,  as  the  cell  size  increases 
to  become  much  larger  than  the  mean-free-path  of 
the  gas. 


3.3  Under-Resolved  Computations 

Given  the  high  over-collision  ratio  in  a  continuum 
application  of  the  DSMC,  one  would  expect,  as  has 
been  argued  above,  that  the  particles  in  each  cell 
will  approach  a  Maxwellian  distribution.  This  is 
confirmed  in  figure  6  which  shows  the  velocity  dis¬ 
tribution  from  an  over-collided  case  {Kc  =  0.0095). 
For  this  case,  the  distribution  was  sampled  at  the 
same  grid  point  and  same  time-level  as  before  (fig¬ 
ure  4)  which,  with  the  current  scaling  corresponds  to 
y  =  245,  t  =  4000.  In  contrast  to  figure  4  however, 
we  see  that  there  is  no  well-defined  structure  to  the 
perturbations  and  that  the  magnitude  of  the  devia¬ 
tions  from  Maxwellian  are  much  smaller.  This  result 
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Figure  6:  Distribution  of  velocity  perturbations  for 
an  over-collided  (under-resolved)  computation,  y  = 
245,  t  =  4000,  The  reference  (Maxwellian)  distri¬ 
bution  has  been  subtracted  and  the  four  quadrants 
have  been  collapsed  to  one  utilizing  the  known  anti¬ 
symmetry  of  the  Chapman-Enskog  distribution. 


Figure  7:  DSMC-computed  viscosity  versus  cell 
Cell  size:  1/Kc  =  Ax.  Under  the  current  non- 
dimensionalization,  the  correct  value  of  is  0.64. 


is  one  we  had  predicted  for  the  over-collided  DSMC 
computations,  but  it  contains  an  implication  that 
is  somewhat  more  subtle.  If  each  cell  in  the  over¬ 
collided  computation  contains  an  equilibrium  distri¬ 
bution,  then  we  are  no  longer  solving  the  Navier- 
Stokes  equations  (which  correspond  to  a  weak  per¬ 
turbation  from  the  Maxwellian  state  of  the  gas)  but 
rather  the  Euler  equations  which  correctly  represent 
an  ideal  gas  in  a  perpetual  state  of  equilibrium.  In¬ 
deed,  the  closer  we  drive  each  cell  to  equilibrium, 
the  more  “ideal”  the  fluid  should  become.  If  this 
is  in  fact  true,  then  the  computed  viscosity  of  the 
Rayleigh  layer  should  go  to  zero  as  the  cell  Knudsen 
number  decreases  and  the  effect  of  the  accommoda¬ 
tion  at  the  solid  surface  should  become  confined  to 
a  thinner  and  thinner  layer,  asymptoting  towards  a 
vorticity  sheet  singularity  at  y  =  0. 

However,  instead  of  a  decrease  in  viscosity  as 
Kc  — ^  0,  we  actually  observe  the  opposite  and,  as 
figure  7  indicates,  the  viscosity  of  the  gas  computed 
from  the  DSMC  computations  increases  as  a  func¬ 
tion  of  the  cell  size.  Ax  =  l/RTc-  For  these  compu¬ 
tations,  the  time  step  was  increased  in  proportion  to 
the  cell  size  so  that  the  CFL  number  was  maintained 
at  0.33  (selected  computations  were  performed  with 
varying  CFL  numbers  yielding  almost  identical  re¬ 
sults.  Three  sets  of  data  are  actually  reported  in 
figure  7.  In  the  first  set  (denoted  by  circles),  the 
number  of  collisions  computed  in  the  collision  phase 
of  each  time  step  was  not  limited.  The  second  series 
of  results,  denoted  by  stars,  a  collision  limiter  of  5 
was  enforced.  Finally,  the  data  indicated  by  crosses 
represents  computations  which  employed  a  collision 
limiter  of  1.  (A  collision  limiter  value  of  one  implies 
that  Np  —  Nc  regardless  of  the  value  of  Np  called  for 
by  the  NTC  equation  (4).  At  low  values  of  l/Kc,  we 
see  that  the  correct  value  of  z/  is  computed.  How¬ 
ever,  as  the  cell  size  rises  above  about  10,  the  value 
of  the  viscosity  begins  to  rise.  The  computed  value 
of  viscosity  appears  to  rise  linearly  over  a  very  wide 
range  of  Kc> 

Upon  reflection,  the  reason  for  this  increase  in  vis¬ 
cosity  is  clear  and  stems  from  the  nature  of  the  par¬ 
ticle  simulation.  At  low  Mach  numbers,  the  random 
thermal  velocity  of  the  gas  will  always  result  in  par¬ 
ticles  moving  up  and  down  in  the  y-direction  even 
though  there  is  no  net  motion  normal  to  the  sur¬ 
face.  Since  the  DSMC  chooses  collision  partners  in 
a  given  cell  without  regard  to  their  location,  momen¬ 
tum  in  the  cell  is  uniformly  diffused  to  neighboring 
cells  within  a  few  time-steps.  In  a  properly-resolved 
computation,  where  the  cell  dimension  is  approxi¬ 
mately  one  mean-free-path,  this  diffusion  is  molec- 
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Figure  8:  Schematic  of  Free-Shear  flow  computa¬ 
tional  geometry. 

ularly  “correct”  and  results  in  the  correct  physical 
viscosity  of  the  fluid  (helped,  of  course,  by  the  appro¬ 
priate  choice  of  VHS  collision  exponent).  However, 
in  a  under-resolved  computation,  the  disregard  for 
a  particle’s  location  in  a  cell  results  in  the  dijffusion 
of  momentum  due  to  thermal  motion  that  is  far  in 
excess  of  the  physical  viscosity  and  it  is  this  arti- 
flcial  viscosity  that  is  exhibited  in  these  results.  In 
many  cases,  this  thermal  motion  would  not  present  a 
problem.  However  in  the  presence  of  a  strong  mean 
velocity  gradient  (as  in  the  case  of  the  wall-driven 
Rayleigh  problem),  the  thermal  motion,  coupled  to 
the  mean  shear,  results  in  an  artificial  Newtonian 
viscosity,  inversely  proportional  to  the  cell  Knudsen 
number,  Kc. 

This  argument  suggests  that  the  artificial  viscos¬ 
ity  is  most  apparent  at  low  Mach  number,  and  that 
as  the  Mach  number  increases,  this  phenomenon  will 
abate  since  the  relative  importance  of  thermal  fluc¬ 
tuations  decrease.  For  the  Rayleigh  problem  this 
was  not  found  to  be  the  case,  but  this  was  ex¬ 
plained  by  the  restriction  of  one- dimensional  flow 
which  was  enforced  in  this  computation  by  wrap¬ 
around  boundary  conditions  which  simply  recycle 
particles  back  into  the  domain  where  they  can  con¬ 
tinue  to  diffuse  in  the  y-direction.  A  more  careful 
assessment  of  Mach  number  effects  is  given  in  the 
following  section. 

3.4  Boundary-Free  Shear  Flows 

In  order  to  verify  this  argument,  it  is  necessary 
to  consider  a  more  complex  geometry  and  we  show 
results  from  a  truly  spatially  developing  flow  -  a 
boundary-free  shear  flow,  illustrated  in  figure  8. 
Here  two  streams  at  different  velocities  are  intro¬ 
duced  at  a:  =  0.  the  upper  stream  is  moving  at 
Ma  =  2.2  while  the  lower  stream  is  moving  at 
Ma  =  1.6. 


Figure  9;  Contours  of  stream  wise  velocity  in  a  well- 
resolved  DSMC  computation  of  a  supersonic  mixing 
layer.  The  x  and  y  dimensions  are  normalized  by  A. 

Figure  9  shows  contours  of  velocity  in  the  x  —  y 
plane  for  a  well-resolved  calculation.  Note  that  the 
singularity  at  a:  =  0  quickly  diffuses  and  a  slowly 
thickening  shear  layer  results.  By  integrating  the 
stream  wise  velocity  across  the  layer,  we  can  charac¬ 
terize  the  flow  by  the  momentum  thickness: 

9  =  -  U){U  -  U2)dy  (13) 

where  0  has  been  normalized  by  the  upper  stream 
velocity,  Ui  and  the  computational  domain  height, 

H.  The  momentum  thickness  is  shown  in  figure  10  as 
a  function  of  x  for  several  computations  at  different 
RTc-  The  theoretical  evolution  of  a  laminar  shear 
layer  [6]  has  a  square- root  growth  of  the  momentum 
thickness: 

9  (X  y/x  (14) 

which  is  well  captured  by  all  of  the  test  cases.  How¬ 
ever,  on  comparing  the  baseline  (well-resolved)  case 
with  the  other  two  cases  (scaling  factors  of  20  and 
40  respectively).  We  would  expect  that,  were  the 
computation  to  replicate  the  correct  level  of  molec¬ 
ular  viscosity  in  the  flow,  the  effect  of  a  twenty¬ 
fold  increase  in  the  linear  dimensions  will  result  in  a 
y/^  =  4.47  decrease  in  the  shear  layer  at  the  same 
value  of  x/H.  In  contrast,  we  observe  that  the  shear 
layer  only  decreases  by  a  factor  of  approximately 

I. 3,  and  that  this  value  is  unchanged  when  the  scal¬ 
ing  factor  is  increased  from  20  to  40  (relative  to  the 
baseline  case).  This  “saturation”  is  due  to  the  fact 
that,  for  this  Mach  number,  the  grid  has  imposed 
its  own  level  of  viscosity  and  that  further  increases 
in  the  scaling  factor  have  no  further  effect  on  the 
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Figure  10:  Development  of  mixing  layer  momentum 
thickness,  0,  as  a  function  of  downstream  distance, 
X.  Each  line  represents  a  different  value  of  Kc  and 
free-stream  Mach  number.  The  a;-direction  is  scaled 
by  the  domain  height,  H  which  ranges  from  60  A  to 
2400. 
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Figure  11:  Geometry  and  Euler  solution  (Mach  con¬ 
tours)  for  Ni’s  Bump.  The  bump  is  20%  of  the  duct 
height.  The  Euler  solution  was  computed  on  a  97 
by  30  grid  with  an  inlet  Mach  number  of  3.28. 


Figure  12:  Contours  of  Mach  number  for  a  DSMC 
computation  in  a  125  by  30  duct.  The  bump  dimen¬ 
sions  are  50  by  6  (All  lengths  normalized  by  A).  The 
inlet  Mach  number  is  3.28  and  the  chord  Reynolds 
number  is  243. 


results  (i.e.  the  “physical  viscosity  has  fallen  below 
the  grid-induced  artifical  viscosity). 

When  the  reference  Mach  number,  Ui/is  increased 
(Keeping  the  velocity  ratio  constant),  we  see  that  the 
shear  layer  thins  further  (still  not  reaching  its  cor¬ 
rect  value,  but  closer  nevertheless).  This  improve¬ 
ment  is  due  to  the  lower  relative  importance  of  ther¬ 
mal  (random)  fluctuations  which  result  in  a  lower 
artifical  viscosity.  It  should  be  noted  that,  as  in  the 
Rayleigh  problem  presented  in  the  previous  section, 
the  artificial  viscosity  was  observed  to  be  relatively 
insensitive  to  the  use  of  a  collision  limiter  and  to  the 
choice  of  the  CFL  number.  In  fact,  it  was  observed 
that  a  collision  limiter  did  thin  the  layer  very  slightly. 
This  was  speculated  to  be  due  to  the  fact  that,  with 
a  collision  limiter,  there  is  slightly  less  “uniformity” 
within  each  cell  since  there  are  less  collisions  and  so 
momentum  is  not  being  smeared  out  quite  so  quickly 
across  the  shear  layer. 

4  Euler  Simulations 

The  results  of  the  preceding  section  illustrate  the 
two  trends  that  compete  in  the  DSMC  computa¬ 
tion  as  the  cell  Knudsen  number  decreases.  On 
the  one  hand,  the  high  value  of  the  over-collision 
ratio,  Np/Nc,  drives  the  cell  particles  towards  a 
Maxwellian,  or  equilibrium,  distribution.  If  this 
were  the  only  effect  present,  the  DSMC  method 


would  formally  solve  the  Euler  equations  which  are, 
after  all,  the  governing  dynamic  equations  for  an 
equilibrium  gas.  However,  a  competing  effect  is  that 
the  essential  nature  of  the  DSMC  technique  -  the 
random  motion  of  particles  coupled  to  the  disregard 
for  location  of  two  collision  partners  -  creates  and  ar¬ 
tificial  viscosity  even  if  the  cell  state  is  Maxwellian. 

From  this,  we  might  expect  that  the  DSMC  code 
should  be  optimal  in  modeling  a  high  Mach  num¬ 
ber,  in  viscid  flow  and  fail  most  dramatically  when 
attempting  a  low-Mach  number  viscous  flow.  The 
latter  has  already  been  demonstrated,  but  we  now 
present  results  to  confirm  our  predictions  of  the  ap¬ 
plicability  of  the  DSMC  method  to  high  Mach  num¬ 
ber  Euler  simulations. 

The  test  case  computed  in  this  investigation  is  a 
common  test  case  for  inviscid  CFD  methods  known 
as  “Ni’s  Bump”  [4]  in  which  the  geometry  consists 
of  a  two-dimensional  duct  with  a  circular-arc  bump 
placed  on  the  floor.  For  the  present  work  we  have 
used  a  bump  with  height  equal  to  20%  of  the  duct 
height.  The  geometry  and  a  traditional  Euler  equa¬ 
tion  solution  is  shown  in  figure  11.  Here  the  inlet 
Mach  number  is  3.28.  A  typical  run  for  these  cases 
contains  about  40,000  particles  in  a  grid  of  about 
2,000  cells.  The  solution  takes  about  10  minutes  to 
reach  statistical  equilibrium.  Adequate  sampling  of 
statistics  takes  one  or  two  more  hours  of  CPU  time 
(on  an  SGI  Indigo). 
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Figure  13;  Contours  of  Mach  number  for  a  DSMC 
computation  in  a  125,000  by  30,000  duct.  The  bump 
dimensions  are  50,000  by  6,000  (all  lengths  normal¬ 
ized  by  A),  The  inlet  Mach  number  is  3.28  and  the 
chord  Reynolds  number  is  243000. 

4.1  Results 

Figure  12  shows  contours  of  Mach  number  for  a 
DSMC  computation  of  Ni’s  bump  in  which  the  duct 
geometry  is  125  by  30  A.  Although  we  are  using 
specular  walls  and  thus  not  attempting  a  viscous  so¬ 
lution  to  the  problem  it  is  instructive  to  know  that 
the  effective  Reynolds  number  of  the  DSMC  com¬ 
putation  (based  on  the  inlet  Mach  number  and  the 
bump  length)  is  243.  In  comparing  the  solution  to 
the  Euler  flow  computation  one  sees  that  the  ba¬ 
sic  solution  is  faithfully  reproduced  with  the  excep¬ 
tion  that  the  shock  thickness  in  the  DSMC  result  is 
overly  large.  This  is  to  be  expected  since,  although 
there  are  no  viscous  effects  due  to  the  boundaries, 
the  DSMC  method  is  accurately  solving  the  shock 
structure,  where  viscous  effects  are  strong  and,  at 
this  low  Reynolds  number,  the  shocks  are  accord¬ 
ingly  thick.  Despite  this,  all  of  the  features  of  the 
flow  are  accurately  predicted  including  the  shock  an¬ 
gles,  the  leading  and  trailing  edge  shocks,  the  shock 
reflections  and  the  shock-shock  interaction  towards 
the  rear  of  the  duct. 

Figure  13  shows  the  same  computation,  now 
scaled  by  a  factor  of  1,000  so  that  the  dimensions  of 
the  duct  are  now  125,000  by  30,000  (approximately 
9  by  2  millimeters  at  atmospheric  conditions).  For 
this  case,  a  collision  limiter  of  one  was  used.  The 
grid  resolution  was  also  doubled  in  each  direction 
(compared  to  the  last  case)  to  approximately  8,000 
cells.  This  was  necessary  since  the  shock  thickness 
was  observed  to  reach  a  minimum  value  regardless 
of  the  cell  Knudsen  number.  As  with  the  shear  layer 
case  presented  above,  this  problem  was  identified  as 
being  an  issue  of  grid  resolution  and  it  was  found 
that,  in  common  with  standard  CFD  practice,  the 
computed  flow  converged  to  a  grid-independent  solu¬ 
tion  through  successive  refinements  in  the  grid  (the 
number  of  test  particles  was  also  increased  although 
the  average  number  of  particles  per  cell  was  smaller 
than  in  previous  case). 


When  compared  with  the  previous  case,  the  shock 
thickness  is  much  smaller,  although  still  larger  than 
“physical”.  This  artificially  thick  shock  is  due  to 
grid  resolution  and  the  inherent  artificial  viscosity 
inherent  in  the  continuum  application  of  the  DSMC 
solution.  However,  this  is  not  to  say  that  the  so¬ 
lution  is  invalid.  It  is  common  practice  in  contin¬ 
uum  CFD  methods  for  solving  the  Euler  Equations 
to  explicitly  introduce  an  artificial  viscosity  (usually 
second-order)  in  order  to  stabilize  the  solution  near 
steep  gradients  such  as  shocks.  The  effect  of  the 
artificial  viscosity  is,  as  in  the  DSMC,  to  thicken 
shocks.  The  DSMC  does  not  require  any  explicit  in¬ 
troduction  of  artificial  viscosity,  it  provides  its  own 
as  demonstrated  by  these  results.  Additionally,  no 
oscillations  in  the  solution  are  observed  in  the  re¬ 
gions  surrounding  the  shock. 

5  Conclusions 

The  discussion  and  results  presented  here  demon¬ 
strate  the  issues  associated  with  scaling  of  a  DSMC 
method  to  continuum  flows.  Two  effects  are  ob¬ 
served  as  the  cell  Knudsen  number  decreases.  The 
first  is,  as  identified  by  Bartel  et  al  [1],  that  each 
cell’s  particles  approach  an  equilibrium  state.  The 
implication  of  this,  and  something  not  previously  re¬ 
alized,  is  that  that  the  DSMC  technique  is  formally 
solving  the  Euler  equations,  not  the  Navier-Stokes 
equations  in  these  regions.  Computationally,  the 
continuum  problem  is  made  tractable  by  employ¬ 
ing  a  collision  limiter  sufficiently  large  so  that  the 
cell  does  approach  equilibrium,  but  small  enough  to 
retain  computational  efficiency. 

A  second  implication  of  the  use  of  DSMC  tech¬ 
niques  for  continuum  methods  is  that  the  thermal 
fluctuations  of  the  particles  in  the  computation,  cou¬ 
pled  to  the  physically  large  cell  size  result  in  an 
artificial  viscosity  which  is  proportional  to  the  cell 
Knudsen  number  and  decreases  with  the  Mach  num¬ 
ber  (as  the  thermal  fluctuations  become  less  impor¬ 
tant).  This  is  a  result  of  the  DSMC  collision  algo¬ 
rithm  which  tends  to  distribute  momentum  and  en¬ 
ergy  uniformly  over  the  entire  cell,  particularly  when 
the  number  of  collision  pairs  during  each  time  step  is 
large  (because  of  the  value  of  Kc) .  This  artifical  vis¬ 
cosity  can  be  regarded  as  dangerous  if  one  is  trying 
to  solve  the  Navier-Stokes  equations  at  low  Knudsen 
numbers,  but  may  be  regarded  as  an  asset  for  solving 
the  Euler  equations  since  it  acts  as  a  built-in  self- 
adjusting  artificial  viscosity  for  broadening  poorly 
resolved  flow  features  such  as  shock  structures. 
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Given  the  advanced  nature  of  continuum  CFD, 
the  motivation  for  using  a  DSMC  method  for  solv> 
ing  a  continuum  problem  might  seem  weak.  How¬ 
ever,  to  our  minds  there  are  several  reasons  in 
which  this  might  be  attractive,  including:  (i)  mixed 
regime  flows  in  which  both  low,  moderate  and  even 
high  Knudsen  regimes  exist  simultaneously,  and  (ii) 
problems  for  which  a  particulate  description  is  ad¬ 
vantageous  -  particle  transport,  multi-species  flows, 
etc.  For  these  applications,  DSMC  approaches  of¬ 
fer  many  advantages.  Lastly,  their  natural  advan¬ 
tages  on  massively  parallel  architectures  might  make 
DSMC  techniques  competitive  with  traditional  CFD 
for  some  applications  and  computational  implemen¬ 
tations. 

Several  techniques  might  be  appropriate  for  accu¬ 
rate  computations  of  continuum  flows  using  particle 
techniques  such  as  the  DSMC,  One  suggestion  is  to 
recognize  that  the  fluid  is  well  described  by  a  well- 
known  distribution  function  (either  a  Maxwellian  or 
a  Chapman-Enskog  distribution)  and  then,  rather 
than  go  through  collision  routines,  simply  assign  the 
particles  in  each  cell  by  sampling  from  a  distribution 
obtained  from  that  cell.  This  “bootstrap”  method 
would  avoid  the  collision  phase  of  the  DSMC  and  re¬ 
place  it  by  a  sampling  phase.  At  this  point  it  is  not 
clear  if  this  would  be  computationally  cheaper  and 
it  is  additionally  not  clear  that  it  would  solve  the  ar- 
tifical  dissipation  due  to  thermal  fluctuations.  The 
idea  is  merely  presented  as  a  possibility  for  contin¬ 
uum  adaptations  to  the  DSMC  approach  warranting 
further  investigation. 
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Nomenclature 


Symbol 

Description 

Section 

a 

speed  of  sound 

2.7.2 

C 

most  probable  thermal  velocity 

2.1 

d 

mean  thermal  velocity 

2.8.1 

Cf 

relative  speed  of  colliding  partners 

1.4.3 

F 

tangential  momentum  accommodation  coeff.  (l=full  accom.) 

2.7.1 

H 

channel  height 

3.1 

Ht 

nozzle  throat  height 

3.3 

k 

Boltzmann  constant,  k  =  1.380658  x  10“^^  J/K 

2.2 

Kn 

Knudsen  number 

1.2 

Kric 

cell  Knudsen  number 

4.1 

Kuo 

outlet  Knudsen  number 

3.1 

L 

channel  length 

3.1 

C 

scale  length  of  flow  gradients 

1.2 

m 

molecular  mass 

2.2 

rrir 

reduced  mass  of  colliding  pair,  =  mim2/{mi  +  m2) 

2.2 

Ma 

Mach  number 

4.2.1 

n 

number  density 

1.4.3 

n 

reference  simulation  number  density 

4.1 

N 

total  number  of  particles 

1.4.2 

Nc 

number  of  particles  in  a  cell 

1.4.3 

Np 

number  of  collision  pairs  to  be  sampled 

2.1 

P 

pressure 

2.7.2 

3 


V 

normalized  pressure,  V  =  plPoutiet 

3,1 

R 

gas  constant,  R  =  k/m 

4.2.1 

Rf 

random  fraction,  0  <  Rj  <  1 

2.7.2 

t 

time 

4.2.1 

T 

temperature 

2.2 

'U,  V,  w 

macroscopic  (mean)  velocity  components 

2.7.2 

u',  v',  w' 

microscopic  (thermal)  velocity  components 

2.7.2 

macroscopic  velocity  component  normal  to  I/O  face 

2.7.2 

Us 

similarity  speed 

3.1 

Uo 

wall  speed  in  Rayleigh  problem 

4.2.1 

Uu  U2 

shear  layer  stream  velocities 

4.2.2 

Vc 

cell  volume 

4.1 

x.y 

spatial  coordinates 

3.1 

1 

specific  heat  ratio 

2.7.2 

At 

time  step 

4.1 

Ax 

typical  cell  dimension 

4.1 

V 

exponent  in  the  inverse  power  law  model 

2.2 

e 

momentum  thickness 

4.2.2 

A 

mean  free  path 

1.2 

y 

coefficient  of  viscosity 

3.1 

V 

kinematic  viscosity  coefficient 

4.2.1 

Vp 

particulate  collision  rate 

2.8.1 

similarity  variable  for  Rayleigh  flow,  ^ 

4.2.1 

a 

cross  section 

1.4.3 

U) 

exponent  in  the  variable  hard  sphere  model 

2.2 
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Chapter  1 


Introduction 


The  direct  simulation  Monte  Carlo  method  (DSMC)  is  a  particle-based  numerical 
fluid  modeling  technique  pioneered  by  G.  A.  Bird.[l]  This  particulate  nature  enables 
it  to  remain  valid  where  traditional,  continuum-based,  computational  fluid  dynamics 
(CFD)  techniques  break  down  due  to  flow  rarefaction.  To  date,  most  work  using 
this  method  has  centered  on  problems  related  to  vehicles  operating  at  high-altitudes, 
where  continuum-based  methods  become  inaccurate  due  to  the  low  ambient  density. 
With  the  advent  of  micro-electro-mechanical  systems  (MEMS),  however,  flows  around 
devices  with  micron-scale  features  have  also  become  important.  Because  their  sharp 
gradients  often  cause  continuum  methods  to  fail,  even  at  standard  conditions,  DSMC 
is  an  attractive  tool  for  these  flows  as  well. 

This  report  explores  DSMC’s  application  to  geometries  related  to  MEMS  devices. 
With  this  task  in  mind,  the  current  chapter  presents  an  introduction  to  MEMS,  non¬ 
continuum  flows,  particle  methods,  and  DSMC.  The  algorithm  developed  for  this 
work  is  then  outlined  in  Chaper  2.  Non-continuum  flows,  the  traditional  realm  of  this 
method,  are  treated  in  Chapter  3  to  evaluate  theoretical  models  as  well  as  explore 
regions  for  which  no  tractable  models  exist.  Chapter  4  then  examines  scaling  issues 
inherent  to  efficiently  applying  DSMC  to  general  continuum  flows,  with  the  intention 
of  expanding  its  useful  range.  Finally,  Chapter  5  presents  conclusions  which  may  be 
drawn  from  this  work  and  points  to  several  features  which  remain  to  be  explored  in 
this  area. 


1.1  MEMS 

Micro-electro-mechanical  systems  (MEMS),  are  very  small  (micron-scale)  sensors  and 
actuators  manufactured  with  techniques  similar  to  those  used  for  integrated  circuit 
(IC)  chips.  They  are  the  subject  of  increasingly  active  research  in  a  widening  field 
of  disciplines.  Applications  for  these  devices  range  from  consumer  products  (airbag 
triggers,  micro-mirror  displays),  to  industrial  and  medical  tools  (microvalves,  micro¬ 
motors),  to  instrumentation  (micro  pressure  sensors,  micro  shear-stress  sensors).  [2] 
MEMS  have  many  advantages  over  their  macro-scale  counterparts,  where  such 
counterparts  even  exist.  First,  because  these  devices  are  fabricated  in  a  manner 
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similar  to  IC  chips,  they  are  extremely  inexpensive  to  manufacture  in  large  quantities. 
Second,  the  technology  for  such  production  is  quite  mature.  Very  precise  specification 
of  the  geometry,  far  beyond  that  possible  with  macro-scale  fabrication  techniques,  is 
therefore  routine  and  a  high  degree  of  control  over  material  properties  is  available. 
Third,  their  small  size  and  mass  make  them  attractive  where  space  is  at  a  premium 
or  weight  is  limited.  Finally,  their  minimal  inertia  allows  them  to  react  very  quickly, 
enabling  the  creation  of  actuators  and  sensors  with  frequency  responses  previously 
unthinkable  for  mechanical  systems. 

The  small  size  of  MEMS  poses  unique  challenges  in  the  design  phase,  however. 
While  the  mechanical  properties  of  micromachined  materials  are  reasonably  well- 
studied,  fluid  effects  at  micron  scales  are  not.  These  effects,  such  as  film  damping  of 
resonant  structures,  heat  transfer  in  mass  flow  sensors,  and  unsteady  pressure  fields 
around  microvalves,  for  example,  must  be  understood  if  the  full  potential  of  these 
devices  is  to  be  realized. 


1.2  Non- Continuum  Flows 

The  vast  majority  of  computational  and  analytical  tools  for  studying  fluid  behavior 
are  based  on  the  Euler  or  Navier-Stokes  equations.  An  important  underlying  assump¬ 
tion  of  these  equations  is  that  the  fluid  may  be  treated  as  a  continuum,  rather  than  as 
a  collection  of  discrete  particles,  as  is  done  in  the,  more  difficult,  Boltzmann  equation. 
This  allows  the  transport  terms  to  be  calculated  using  macroscopic  variables,  such  as 
temperature,  rather  than  microscopic  variables,  such  as  molecular  velocity  distribu¬ 
tion  function,  yielding  an  expression  which  is  much  more  amenable  to  solution,  both 
analytically  and  numerically.  Unfortunately,  this  approximation  becomes  inaccurate 
as  the  scale  length  of  the  flow  gradients  (E)  approaches  the  average  distance  trav¬ 
eled  by  a  particle  between  collisions  (the  mean  free  path,  A),  which  occurs  for  many 
MEMS-related  flows.  The  ratio  of  these  quantities  is  known  as  the  Knudsen  num¬ 
ber  {Kn  =  \/C)  and  is  used  to  indicate  the  degree  of  flow  rarefaction.  The  Navier- 
Stokes  equations  neglect  rarefaction  effects  and  are  therefore  only  strictly  accurate 
for  vanishingly-small  Kn. 

As  is  the  case  for  other  non-dimensional  numbers  in  fluid  dynamics,  such  as  the 
Mach  and  Reynolds  numbers,  the  type  of  analysis  appropriate  for  a  particular  flow 
is  dictated  by  its  Knudsen  number.  Consequently,  the  Kn  domain  (0  <  Kn  <  oo) 
is  often  divided  into  four  flow  regimes. [3]  For  Kn  <  0.01,  known  as  the  ‘contin¬ 
uum’  regime,  the  Navier-Stokes  equations,  as  commonly  expressed,  are  applied.  For 
0.01  <  Kn  <  0.1,  known  as  the  ‘slip-flow’  regime,  the  Navier-Stokes  equations  are 
applied  with  the  usual  no-slip  wall  boundary  condition  replaced  by  a  slip-flow  condi¬ 
tion  (detailed  in  Section  3.1).  For  0.1  <  Kn  <  3,  known  as  the  ‘transition’  regime, 
the  flow  is  too  rarefied  for  Navier-Stokes-based  analysis,  but  not  rarefied  enough  to 
apply  the  collisionless  Boltzmann  equation.  The  full  Boltzmann  equation  is  therefore 
prescribed.  For  Kn  >  3,  known  as  the  ‘free  molecular’  regime,  the  flow  is  sufficiently 
rarefied  to  allow  molecular  collisions  to  be  neglected.  The  collisionless  Boltzmann 
equation  is  therefore  applied. 
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For  many  MEMS.  Kn  is  driven  from  the  continuum  regime  by  their  extremely 
small  feature  size,  which  is  often  comparable  to  A,  even  at  standard  conditions.  The 
gap  between  the  sensing  plate  and  the  substrate  on  a  floating  element  shear  stress 
sensor,  for  example,  is  typically  1-2  microns.  [4]  The  mean  free  path  of  air  at  standard 
conditions  is  approximately  60  nm.  This  places  Kn  in  the  slip  flow  regime,  or  even 
the  transition  regime  for  lighter  gases  or  different  conditions.  Non-continuum  effects 
in  the  gap,  neglected  in  traditional  analyses,  may  therefore  have  a  significant  impact 
on  sensor  operation.  As  a  result,  new  models  and  techniques  must  be  developed  to 
correctly  describe  the  behavior  of  the  fluid  in  and  around  these  devices. 

1.3  Particle  Methods 

Particle  methods,  such  as  molecular  dynamics  (MD),  particle-in-cell  (PIC),  and  direct 
simulation  Monte  Carlo  (DSMC)  are  attractive  tools  for  the  study  of  rarefied  gas- flows 
because  they  lack  continuum  assumptions.  These  techniques  model  gas  behavior  by 
tracking  the  interaction  of  computational  particles,  each  with  a  position,  a  velocity, 
an  internal  energy,  etc.,  mimicking  the  discrete  molecular  nature  of  the  actual  flow. 
This  strategy  differs  considerably  from  that  of  traditional  CFD,  which  numerically 
solves  differential  field  equations  formulated  to  describe  fluid  behavior  in  terms  of 
macroscopic  variables. 

Particle  methods  are  intuitively  attractive  because  fully  physical  simulations  would 
be  devoid  of  assumptions  and  would  therefore  be  valid  for  all  flow  regimes  and  geome¬ 
tries.  Unfortunately,  a  fully-physical  simulation,  even  for  a  simple  problem,  typically 
requires  computational  power  several  orders  of  magnitude  greater  than  is  currently 
available.  In  addition,  it  is  arguable  that  a  complete  understanding  of  intermolecular 
and  molecule-surface  interactions  is  not  yet  reached,  so  some  detail  would  necessarily 
be  neglected  if  such  a  calculation  were  even  attempted. 

Fortunately,  many  features  of  molecular  behavior  have  negligible  influence  on  most 
engineering  problems.  The  key  to  an  efficient  simulation  is  therefore  to  include  the 
minimum  level  of  complexity  required  to  correctly  reproduce  the  important  features 
of  a  given  flow.  Consequently,  all  current  particle  methods  make  simplifications  in 
the  quantity,  movement,  and/or  interactions  of  their  computational  particles.  The 
form  of  these  simplifications  differentiates  between  the  various  techniques. 

1.4  DSMC 

DSMC  is,  by  far,  the  most  popular  particle  method  for  the  analysis  of  collisional 
flows,  ie.  flows  for  which  intermolecular  collisions  significantly  affect  fluid  behavior. 
It  is  called  a  simulation  (rather  than  a  solution)  scheme  because  it  was  originally 
formulated  to  capture  the  important  physical  features  of  the  flow,  not  to  solve  a  par¬ 
ticular  set  of  equations.  Nevertheless,  Nanbu  later  showed  that  the  various  techniques 
solve  either  the  Kac  Master  equation  or  the  Boltzmann  equation. [5]  As  a  result,  some 
current  algorithms  were  subsequently  derived  from  these  equations,  rather  than  from 
physical  arguments. 
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DSMC  has  undergone  considerable  development  by  many  researchers  for  more 
than  two  decades.  Consequently,  a  considerable  number  of  variations  from  the  original 
algorithm  now  exist.  A  number  of  features  nonetheless  remain  common  to  most 
implementations.  These  defining  features  are  associated  with  the  simplifications  made 
to  the  physical  situation  and  differentiate  DSMC  from  other  particle  methods. 

1.4.1  Particle  Quantity 

The  number  of  particles  in  the  simulation  must  be  reduced  due  to  memory  constraints 
as  well  as  CPU  time  considerations.  For  example,  a  cubic  centimeter  of  gas  at  stan¬ 
dard  conditions  contains  approximately  2.69  x  10^^  molecules.  To  represent  only 
position  and  velocity  in  3-dimensional  space  for  these  particles  with  single- precision, 
floating  point  numbers  would  require  a  total  memory  of  6.5  x  10^“^  Mb,  far  beyond  the 
capacity  of  even  modern  supercomputers. 

The  crucial  assumption  made  by  DSMC  developers  in  the  face  of  this  problem 
is  that,  at  any  given  time,  a  number  of  molecules  are  in  virtually  indistinguishable 
microscopic  states.  These  molecules  may  therefore  be  collectively  represented  as  a 
single  ‘computational  particle’,  which  is  then  a  statistical,  rather  than  a  physical, 
entity  because  it  represents  a  small  volume  in  phase  space,  rather  than  an  actual 
molecule. 

The  ratio  of  real  to  computational  particles  is  known  as  the  ‘weight  factor’  and 
may  be  constant  for  a  given  computation  or  vary  with  position  and/or  time.  Variable 
weight  factors  are  useful  in  rapidly  expanding  and  axisymmetric  flows,  where  large 
changes  in  number  density  or  cell  volume  occur  across  the  domain,  to  hold  the  num¬ 
ber  of  particles  in  each  cell  at  a  computationally  efficient  level.  [6]  These  weighting 
schemes  must  be  applied  with  caution  because  particles  are  usually  ‘cloned’  when  they 
move  into  a  region  with  a  smaller  weight  factor.  This  is  problematic  because  cloned 
particles  are  not  statistically  independent.  Introducing  significant  numbers  of  them 
will  therefore  cause  a  larger  scatter  to  be  observed  in  the  results  than  is  expected  for 
the  given  number  of  particles.  In  addition,  their  lack  of  relative  speed  may  cause  a 
non-physical  reduction  of  the  local  temperature. 

The  accuracy  of  the  particle  quantity  approximation  increases  with  the  number  of 
computational  particles  because  each  particle  is  required  to  represent  a  smaller  region 
of  phase  space.  In  addition,  the  statistical  scatter  in  a  given  sample  of  the  flowfield 
is  reduced  through  increased  particle  quantity.  These  factors  are  balanced  by  the 
computational  cost  and  machine  requirements  of  the  simulation,  which  also  increase 
with  particle  quantity.  This  issue  is  treated  in  detail  by  Chen  and  Boyd  in  Ref.  [7]. 

It  is  notable  that  MD  also  uses  a  reduced  number  of  computational  particles  but 
continues  to  treat  them  as  physical  entities.  This  is  justified  by  noting  that  Kn  is 
maintained  (for  flow  similarity  between  the  real  and  simulated  cases)  if  the  particle 
diameter  is  increased  to  hold  the  mean  free  path  constant  as  the  number  density 
is  reduced.  [8]  A  large  number  of  small  particles  may  therefore  be  replaced  with  a 
small  number  of  large  particles  without  affecting  the  important  flow  characteristics. 
This  argument  fails,  however,  when  the  resulting  particle  diameters  become  too  large 
compared  to  the  geometric  feature  size  or  the  inter-particle  spacing. 
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1.4.2  Particle  Movement 


In  a  real  gas,  molecules  move  along  their  current  trajectories  until  they  strike  another 
molecule  or  a  boundary.  This  is  the  procedure  used  for  movement  in  MD,  consis¬ 
tent  with  its  treatment  of  computational  particles  as  physical  entities.  Unfortunately, 
each  particle’s  next  collision  depends  on  its  trajectory  and  position  with  respect  to 
all  other  particles  in  the  domain,  as  well  as  the  boundaries.  In  the  worst  case,  the 
computational  work  of  this  combined  move/collide  operation  therefore  scales  as 
where  N  is  the  total  number  of  computational  particles.  While  clever  ways  of  mini¬ 
mizing  this  work  have  been  proposed,  such  as  storing  the  time  of  next  interaction  for 
all  particle  pairs  and  recalculating  only  when  it  changes  through  the  collision  of  one 
or  both  members[8],  the  computation  quickly  becomes  unmanageable  as  the  number 
of  particles  increases. 

One  of  DSMC’s  defining  assumptions  is  made  in  response  to  this  problem:  if  par¬ 
ticles  are  only  allowed  to  move  for  a  short  time  (some  fraction  of  a  molecule’s  mean 
time  between  collisions),  then  the  movement  and  intermolecular  collision  steps  may 
be  decoupled.  The  calculation  is  therefore  divided  into  ‘time  steps’,  each  consisting  of 
independent  move  and  collide  phases.  Interactions  between  particles  and  the  domain 
boundaries  are  still  handled  as  they  occur  in  the  move  phase,  but  all  intermolecu¬ 
lar  collisions  are  performed  in  the  collide  phase.  This  reduces  movement  from  an 
order  N'^  to  approximately  an  order  N  calculation.  It  may  be  noted  that  this  also 
makes  it  possible  for  two  molecules  to  simultaneously  occupy  the  same  position.  This 
seemingly  troublesome,  non-physical  event  is  not  significant,  however,  because  the 
computational  particles  each  represent  a  range  of  positions  due  to  their  statistical 
nature.  Their  exact  position  on  the  grid  is  therefore  not  important. 


1.4.3  Particle  Collisions 

In  MD,  interparticle  collisions  are  physical  events,  with  the  collision  partners  and 
post-collision  velocities  determined  according  to  the  pre-collision  particle  trajectories 
and  diameters.  These  diameters,  as  mentioned  above,  are  scaled  to  maintain  the 
proper  collision  rate  for  the  reduced  number  of  particles  in  the  calculation. 

For  DSMC,  however,  the  simplifications  made  above  have  eliminated  the  physical 
means  for  determining  these  collisions.  An  alternate  approach  must  therefore  be 
supplied.  This  involves  both  calculating  the  correct  number  of  collisions  to  perform 
and  choosing  proper  partners  for  each  of  them. 


Quantity  Calculation 

Physical  arguments  from  classical  kinetic  theory  may  be  used  to  develop  an  expression 
for  the  number  of  collisions  per  unit  time  per  unit  volume  of  a  gas,  Nc- 


Nc  =  -n^  crCr 
2 


(1.1) 
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where  n  is  the  number  density,  a  is  the  molecular  cross-section  (single-species  flow), 
Cr  is  the  relative  speed  of  the  colliding  partners,  and  the  overbar  signifies  a  mean 
taken  over  all  possible  partners.  Unfortunately,  the  computational  work  to  evaluate 
this  mean  scales  as  the  total  number  of  particles  squared.  Several  ways  of  eliminating 
this  term  have  therefore  been  proposed,  such  as  Bird’s  Time  Counter  (TC)  and  No 
Time  Counter  (NTC)  schemes  and  Baganoff  and  McDonald’s  method. [9] [10] [11]  The 
current  code  uses  the  NTC  scheme,  where  the  local  maximum  acr  encountered  in  the 
calculation  is  used  in  Eq.  1.1  to  determine  the  number  of  particle  pairs  to  consider 
for  collision.  These  pairs  are  then  accepted  with  a  probability  proportional  to  their 

(7  C-p . 

Partner  Selection 

Choosing  potential  collision  pairs  is  best  done  according  to  some  sort  of  ‘close  physical 
proximity’  criterion.  The  most  computationally  convenient  of  these  is  to  divide  the 
domain  into  a  number  of  ‘cells’,  as  is  commonly  done  for  traditional  CFD  calculations. 
Collision  pairs  are  then  chosen  among  particles  in  the  same  cell.  For  this  to  be  a  valid 
approximation,  these  cells  must  be  ‘small’.  Rigorously,  this  means  that  their  linear 
dimensions  should  be  comparable  to  the  mean  free  path,  A.  The  implications  of 
relaxing  this  requirement  to  the  ‘negligible  gradients  across  the  cell’  constraint  used 
in  continuum  CFD,  are  investigated  in  Chapter  4. 

1.4.4  Result  Reporting 

The  cells  required  by  the  collision  scheme  are  also  used  in  reporting  the  results  of 
a  simulation.  This  is  necessary  because  macroscopic  variables,  such  as  pressure  and 
temperature,  are  typically  the  quantities  of  interest  from  a  computation  while  the 
method  functions  in  terms  of  microscopic  variables,  such  as  individual  particle  posi¬ 
tions  and  velocities.  To  determine  the  former  from  the  latter,  the  state  of  all  particles 
in  some  small  volume  surrounding  the  point  of  interest  must  be  sampled.  In  the 
current  code,  the  particles  in  each  cell  are  used  to  calculate  the  macroscopic  variables 
which  are  reported  at  its  centroid. 

Due  to  the  total  particle  quantity  restrictions,  imposed  by  computational  consid¬ 
erations,  and  the  cell  size  constraints,  imposed  by  collision  and  sampling  concerns, 
there  may  be  as  few  as  twenty  computational  particles  in  a  given  cell.  Calculating 
the  macroscopic  variables  from  a  single  sample  will  therefore  yield  unacceptably  large 
uncertainties.  Several  samples  of  each  cell  are  therefore  necessary.  If  steady-state 
information  is  required,  these  samples  may  be  taken  every  few  time  steps^  until  ac¬ 
ceptable  statistical  convergence  is  reached.  If  time-accurate  information  is  required, 
ensemble  averaging  is  performed,  were  the  entire  flow  evolution  is  calculated  repeat¬ 
edly,  sampling  at  the  same  temporal  locations  in  each  iteration. 

A  flowchart  of  a  typical  DSMC  calculation  is  presented  in  Figure  1-1. 


The  samples  would  not  be  statistically  independent  if  they  were  taken  every  time  step 


10 


Figure  1-1:  Flowchart  for  a  typical  DSMC  calculation. 
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Chapter  2 
Algorithm 


Due  to  the  wide  variety  of  geometries  and  flow  conditions  present  in  MEMS  devices, 
flexibility  was  a  primary  goal  in  constructing  the  algorithm  for  this  investigation. 
DSMC  has  an  inherent  advantage  over  traditional  CFD  in  this  regard  because  its 
particulate  nature  makes  it  uniformly  valid  for  all  Kn,  without  regime-specific  modi¬ 
fications.  This  is  especially  important  for  MEMS  work  because  many  devices  contain 
mixed-regime  flows. 

Flexibility  concerns  also  drove  the  particular  implementation  of  DSMC  developed 
for  this  work.  First,  C  was  chosen  for  the  coding  due  to  its  flexible  structure,  dy¬ 
namic  memory  allocation,  and  recursive  function  capability.  In  addition,  the  code  was 
written  to  run  in  non-dimensional  form,  use  unstructured  grids,  and  store  most  data 
locally.  The  details  of  the  resulting  algorithm  are  discussed  in  the  following  sections. 


2.1  Non-Dimensionalization 

All  quantities  in  the  code  are  non-dimensionalized.  This  generalizes  the  construction 
of  calculations  and  facilitates  the  interpretation  of  their  results.  The  normalization 
factors  selected  for  this  purpose  are  relatively  common  in  the  DSMC  community. 

Due  to  the  importance  of  Kn  in  the  subject  geometries,  the  mean  free  path,  at 
some  reference  condition,  is  a  logical  choice  for  non-dimensionalizing  length.  Simi¬ 
larly,  DSMC’s  particulate  nature  points  to  one  of  the  molecular  speeds  (ie.  mean, 
rms,  or  most  probable)  as  the  normalizing  factor  for  velocity.  These  speeds  differ  by 
a  constant  factor  near  unity,  so  the  choice  is  essentially  arbitrary.  The  most  probable 
molecular  speed,  c^,  is  used  in  this  work.  The  quotient  of  the  length  and  speed  nor¬ 
malizations,  A/c'^,  \pK  times  the  mean  collision  time)  is  used  to  non-dimensionalize 
time.  Temperature,  pressure,  and  number  density  are  non-dimensionalized  by  their 
values  at  a  reference  condition. 

To  non-dimensionalize  the  expression  governing  the  number  of  collision  pairs  to 
be  sampled  in  the  NTC  scheme: 


A/p  —  n  (^(7Cr')j 


(2.1) 
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Variable 

Scale  Factor 

Length 

^ref 

Velocity 

d 

TUref 

Time 

^ref 

c! 

^ref 

Temperature 

Tref 

Pressure 

Pref 

Number  Density 

Tire  f 

Collision  Cross-Section 

1 

^ref^ref 

Table  2.1:  Non-Dimensionalization  /  Normalization  Factors 


a  normalizing  factor  must  be  chosen  for  the  collision  cross-section,  a.  Here  it  should  be 
noted  that  the  probability  of  a  particular  molecule  suffering  a  collision  is  proportional 
to  the  product  of  its  cross-section,  its  path  length,  and  the  number  density.  For  a 
non-dimensional  flow  representation  to  be  similar  to  its  dimensional  counterpart,  the 
collision  probabilities  for  comparable  particles  must  match.  As  a  result,  the  scale 
factors  for  collision  cross-section,  length,  and  number  density  must  have  a  product 
of  unity.  The  proper  cross-section  scale  factor  is  therefore  l/(nre/Are/).[12]  This 
conclusion  may  also  be  reached  through  purely  dimensional  arguments  if  the  units 
for  a  are  viewed  as  ‘length-units  squared  per  particle’,  rather  than  simply  ‘length- 
units  squared’. 

The  non-dimensionalization  /  normalization  factors  used  in  the  code  are  summa¬ 
rized  in  Table  2.1. 


2.2  Molecular  Model 

An  important  defining  feature  of  any  DSMC  code  is  its  molecular  model.  This  spec¬ 
ifies  the  collision  cross-section  used  in  Eq.  2.1  as  well  as  the  post-collision  scattering 
law.  There  is,  again,  a  tradeoff  to  be  made  between  physical  realism  and  compu¬ 
tational  efficiency.  The  full  physical  description  of  intermolecular  interaction  is  not 
known  and,  for  most  flows,  not  necessary.  As  a  result,  many  models  have  been  pro¬ 
posed  which  attempt,  with  varying  degrees  of  detail,  to  enable  the  computation  to 
reproduce  the  important  features  of  a  collisional  flow,  such  as  the  viscosity  coefficient 
and  its  temperature  dependence,  without  becoming  unnecessarily  complex.  Exam¬ 
ples  include  the  Inverse  Power  Law  (IPL),  Hard  Sphere  (HS)  and  Maxwell  models  of 
classical  kinetic  theory,  the  Variable  Hard  Sphere  (VHS)  model  of  Bird,  the  Variable 
Soft  Sphere  (VSS)  model  of  Koura  and  Mastumoto,  and  the  Generalized  Hard  Sphere 
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(GHS)  model  of  Hassan  and  Hash.[13][14][15][16] 

The  V'HS  model  is  implemented  in  the  current  code.  This  model  grew  from 
the  classical  HS  representation,  which  assumes  a  constant  collision  cross  section  and 
isotropic  scattering.  For  real  molecules,  however,  the  effective  cross-section  is  reason¬ 
ably  constant  only  at  very  low  temperatures;  at  higher  temperatures,  it  decreases  with 
increasing  translational  kinetic  energy  and  relative  speed  of  the  colliding  partners. [1] 
In  addition,  post-collision  scattering  is  decidedly  non- isotropic.  Bird  noted,  from  nu¬ 
merical  and  analytical  studies,  that  changes  in  collision  cross-section  have  a  strong 
influence  on  the  gas  behavior  while  changes  in  the  scattering  law  do  not.  He  therefore 
constructed  a  model  with  a  Cr-dependent  cross-section,  but  with  isotropic  scattering; 
essentially  creating  hard-sphere  molecules  with  variable  diameters  (hence  the  name). 
This  retains  most  of  the  computational  simplicity  of  the  HS  model,  but  more  accu¬ 
rately  reproduces  the  temperature  dependence  of  the  viscosity  coefficient. 

In  this  model,  an  empirical  constant,  u,  related  to  the  exponent,  rj,  of  the  inverse 
power  law  molecular  force  by 

u;  =  2/(r?-l).  (2.2) 

is  supplied  to  establish  the  working  fluid.  The  collision  cross-section  is  then  given  by 


Cfd  — 


(  \ 
\2(2  —  Uj)kTref^  ) 


(2.3) 


where  k  is  Boltzmann’s  constant  and  is  the  reduced  mass  of  the  colliding  pair, 
TUr  =  mi 7712 /(mi  -I- m2),  which  is  always  ml 2  for  the  single-species  gases  considered  in 
this  work  (m  is  the  molecular  mass).  Subscript  “d”  denotes  a  dimensional  quantity. 

Under  the  normalizations  introduced  in  the  previous  section,  using  the  VHS  mean 
free  path  given  by  Bird: 

Ad  =  {T|Trefr|[^^2{2  -  a;)^^  r(2  -  uj)ndares,i  (2.4) 


the  non-dimensional  collision  cross-section  becomes: 

where  r()  denotes  the  gamma  function. 


(2.5) 


2.3  Grids 

In  response  to  the  stated  goal  of  maximum  flexibility,  the  current  code  was  written 
for  unstructured  grids.  This  enables  it  to  treat  geometries  of  arbitrary  complexity 
without  modification.  In  addition,  several  sophisticated  generation  and  adaptation 
schemes  are  available  for  grids  of  this  type.  Though  the  code  was  written  with  a 
generalized  cell  geometry  in  mind,  all  calculations  discussed  in  Chapters  3  and  4 
use  3-sided  cells  generated  as  a  Delaunay  triangulation  of  points  distributed  in  the 
domain  and  on  the  boundaries.  This  triangulation  was  performed  using  Watson’s 
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algonthm[17],  which  made  possible  the  inclusion  of  ‘point-and-click’  node  placement 
for  pre-calculation  local  refinement. 

A  sample  grid  generated  in  this  manner  is  presented  in  Figure  2-1.  The  geometry 
shown  is  a  channel  with  a  25%  bump,  similar  to  the  case  run  in  Section  4.2.3.  Grid 
refinement  is  demonstrated  at  the  left  edge  of  the  bump. 


Figure  2-1:  Example  of  an  unstructured  grid  with  local  refinement  generated  with 
Watson’s  algorithm. 


2.4  Particle  Movement 

It  is  common  practice  for  DSMC  codes  on  structured  grids  to  displace  all  particles 
a  full  time  step  along  their  trajectories,  then  determine  their  resulting  cell  indices 
through  a  search  or  mathematical  operation.  By  their  nature,  however,  unstructured 
grids  make  this  type  of  scheme  very  difficult.  A  solution  to  this  problem  was  pro¬ 
posed  (though  for  structured  3-D  grids)  by  Dietrich[18]:  perform  particle  movement 
and  current-cell  identification  simultaneously  by  maintaining  knowledge  of  a  parti¬ 
cle’s  current  cell  at  all  points  in  its  trajectory.  To  accomplish  this  computationally, 
a  particle  is  displaced  until  it  contacts  a  face  of  its  current  cell.  It  is  then  passed 
to  an  adjoining  cell,  reflected  from  a  solid  boundary,  or  allowed  to  leave  the  calcula¬ 
tion  through  an  inflow/outflow  edge,  as  specified  by  a  ‘neighbor  identifier’  stored  for 
each  face  of  the  cell.  This  arrangement  is  extraordinarily  flexible  because  the  grid 
need  only  be  composed  of  cells  with  valid  neighbor  identifiers  on  each  of  their  faces. 
These  cells  can,  strictly  speaking,  possess  an  arbitrary  number  of  edges  and  be  of  any 
shape  or  orientation.  From  a  computational  viewpoint,  however,  this  method  is  most 
convenient  if  the  cells  are  required  to  be  convex.  This  allows  the  faces  to  be  treated 
as  infinitely  long  and  the  impacted  face  to  be  found  by  selecting  the  line  that  the 
particle  trajectory  intersects  at  the  earliest  time. 

This  portion  of  the  code  proved  to  be  the  most  sensitive  to  numerical  accuracy 
issues.  Due  to  the  imperfect  machine  resolution  of  position  and  cell  information,  it  is 
difficult  to  calculate  the  proper  destination  cell  for  particles  which  are  close  to  faces  or 
cross  near  nodes.  A  number  of  measures  were  implemented  to  combat  this  problem. 

The  first  of  these  involves  noting  the  face  through  which  a  particle  entered  its  cell. 
This  face  is  then  excluded  from  consideration  as  a  crossing  site  for  the  remainder  of 
the  current  time  step.  Unfortunately,  due  to  the  unstructured  nature  of  the  grid, 
no  mathematical  means  exist  for  identifying  the  entry  face  of  the  new  cell  from  the 
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(known)  exit  face  of  the  old  cell.  A  search  for  a  face  in  the  new  cell  whth  a  neighbor 
identifier  pointing  to  the  particle’s  former  cell  is  therefore  performed. 

From  this  effort,  calculating  a  useless  (and  problematic)  intersection  time  is  avoided. 
This  quantity  is  problematic  because  particles  are  displaced  to  their  intersection  point 
when  passed  to  a  new  cell.  The  intersection  time  in  the  new  cell  for  this  face  should 
therefore  be  zero,  but  will,  in  practice,  be  a  small  number  due  to  machine  resolu¬ 
tion  issues.  If  this  number  is  positive,  this  face  is  likely  to  be  selected  as  the  next 
intersection  because  it  has  the  smallest  time  to  occurrence,  causing  the  particle  to  be 
passed  back  to  its  original  cell,  which  will  return  it,  resulting  in  an  infinite  loop.  The 
computation  and  memory  cost  of  searching  for  and  storing  a  particle’s  entry  face  is 
therefore  justified. 

Another  test  is  included  for  crossings  which  occur  very  close  to  cell  nodes,  such  as 
that  shown  in  Figure  2-2.  In  this  case,  the  particle  will  be  transferred  from  cell  one  to 
two  and  then  to  three  without  difficulty.  A  problem  is  encountered  when  searching  for 
intersections  in  cell  three,  however.  The  face  between  cells  two  and  three  is  rejected 
because  it  was  just  crossed,  but  an  intersection  time  near  zero  is  calculated  for  the 
other  face  which  shares  the  node  skirted  by  the  particle.  This  causes  the  particle  to  be 
(erroneously)  passed  to  cell  four  without  significantly  moving  it  (due  to  the  miniscule 
intersection  time).  In  a  similar  fashion,  cell  four  transfers  the  particle  to  cell  one, 
which  transfers  it  to  cell  two,  and  so  on.  This  results  in  an  unending  cycle  because 
each  of  these  moves  has  a  negligible  duration,  so  the  time  step  is  never  completed. 


Figure  2-2;  Sketch  of  a  particle  crossing  near  a  grid  node 

To  avoid  establishing  this  cycle,  a  ‘direction  test’  is  performed  on  potential  inter¬ 
section  faces.  In  this  test,  the  particle  trajectory  is  projected  onto  the  inward-pointing 
normal  of  the  face  in  question.  If  the  projection  has  a  positive  sign,  then  the  particle 
is  not  moving  in  the  proper  direction  to  strike  this  face  and  the  intersection  is  rejected. 
To  avoid  encumbering  routine  cases  with  this  test,  a  tolerance  is  defined.  The  test 
is  then  only  performed  on  faces  whose  intersection  times  fall  within  this  tolerance  of 
zero. 

A  code  listing  of  the  movement  function  is  provided  in  Appendix  A. 


2.5  Data  Structure 

Many  DSMC  codes  use  a  single  array  to  store  the  data  for  all  particles  in  the  simula¬ 
tion.  A  ‘cross-reference  vector’  is  then  maintained  by  each  cell,  containing  the  indices 
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of  its  particles  in  the  central  array.  Moving  a  particle  from  cell  to  cell  then  consists 
of  simply  transferring  its  index  from  one  cross-reference  vector  to  another.  This  is  a 
very  efficient  arrangement  for  machines  with  rapid  access  to  all  their  memory,  such 
as  supercomputers,  because  very  little  data  must  be  moved  with  the  particle. 

Dietrich  and  Boyd  noted,  however,  that  this  scheme  causes  a  significant  loss  of 
computational  efficiency  on  ‘workstation’  computers.  [19]  These  machines,  unlike  su¬ 
percomputers,  have  a  relatively  slow  main  memory  but  a  very  fast  cache.  Before 
performing  an  operation,  this  cache  is  loaded  with  the  segment  of  main  memory  con¬ 
taining  the  required  data.  If  subsequent  operations  access  only  items  already  in  the 
cache,  they  execute  very  quickly.  Conversely,  if  the  cache  must  be  reloaded  with  new 
data,  a  condition  known  as  a  ‘cache  miss’,  a  considerable  amount  of  time  is  lost  in 
the  process.  The  central  storage  of  particle  data  causes  many  cache  misses  when 
performing  inherently  cell-based  operations,  such  as  collisions  and  sampling,  because 
the  members  of  a  given  cell  are  scattered  through  a  very  large  area  of  main  memory. 

The  alternative  is  clear:  physically  store  a  particle’s  data  in  its  current  cell  instead 
of  in  a  large,  centralized  array.  This  increases  the  expense  of  moving  a  particle  from 
cell  to  cell,  but  greatly  decreases  the  number  of  cache  misses  suffered  by  cell-based 
operations. 

An  algorithm  using  the  trajectory-tracing  particle  movement  scheme  outlined  in 
the  previous  section  is  easily  written  with  entirely  cell-based  operations,  so  consid¬ 
erable  gains  are  possible  from  improved  utilization  of  the  cache.  In  addition,  for 
an  equilibrium  gas  with  a  properly  chosen  cell  size  and  time  step,  less  than  50%  of 
particles  typically  leave  their  cell  during  a  given  time  step,  and,  of  these,  less  than 
10%  typically  leave  their  new  cell.  The  increased  cost  of  inter-cell  movement  in  this 
arrangement  is  therefore  greatly  overshadowed  by  the  increased  speed  of  cell-based 
data  structures. 

This  data  structure  has  the  additional  advantage  of  easy  adaptability  to  a  message¬ 
passing  parallel  environment.  If  a  particle  moves  to  a  cell  on  a  different  processor,  the 
inter-cell  communication  step  is  simply  augmented  with  an  inter-processor  commu¬ 
nication  step.  In  addition,  because  a  cell  is  now  a  complete  data  unit,  containing  its 
geometry,  neighbors,  and  resident  particles,  it  is  easily  passed  as  a  whole  to  another 
processor  as  part  of  a  dynamic  load  balancing  scheme. 

2.6  Inter-Cell  Communication 

Upon  adopting  a  cell-based  data  structure,  an  efficient  means  for  moving  particle 
information  between  cells  is  required.  The  simplest  method,  which  is  similar  to  that 
used  by  Dietrich  and  Boyd  for  communication  between  parallel  domains,  is  to  main¬ 
tain  a  ‘communication  link’  for  each  cell.  This  is  an  array  for  particles  waiting  to 
be  accepted  into  the  cell.  If  a  particle  leaves  its  cell  during  the  movement  step,  it 
is  simply  placed  in  the  communication  link  of  its  destination  cell.  When  a  cell  has 
completed  the  movement  step  for  all  of  its  particles,  it  begins  to  process  its  com¬ 
munication  link:  each  of  the  incoming  particles  is  checked  to  see  if  it  will  remain  in 
the  cell;  if  so,  it  is  moved  to  the  particle  list  for  that  cell,  if  not,  it  is  moved  to  the 
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communication  link  of  its  next  cell.  This  process  continues  until  all  communication 
links  are  empty.  This  arrangement  was  found  to  be  fast,  but  extremely  demanding  of 
memory;  the  communication  links  typically  grew  to  be  about  half  as  large  as  the  par¬ 
ticle  arrays  for  their  cell.  The  total  memory  required  for  the  program  was  therefore 
50%  larger  than  without  links. 

To  circumvent  this  difficulty,  a  new  scheme  was  devised  whereby  particles  are 
simply  marked  with  a  flag  signifying  that  they  are  not  leaving  the  cell,  the  index 
of  their  destination  cell,  or  a  flag  signifying  that  they  have  exited  the  cell  and  their 
position  on  its  particle  list  is  now  vacant.  Once  the  initial  movement  step  is  completed 
for  all  particles  in  all  cells,  a  ‘communication  phase’  is  entered  in  which  particles 
leaving  their  current  cells,  referred  to  as  ‘travelers’,  are  moved  to  their  destination 
cells.  To  accomplish  this,  the  particle  array  in  a  traveler’s  destination  cell  is  checked 
for  empty  positions.  If  one  is  found,  the  particle  is  moved  to  its  new  cell,  marking  its 
former  position  as  vacant.  If  no  vacancies  are  found,  the  destination  cell  is  searched 
for  travelers.  If  a  traveler  is  found,  the  communication  function  is  recursively  called 
to  move  this  particle  to  its  destination  cell,  then  the  original  particle  is  moved  into 
the  space  it  left.  If  neither  vacancies  nor  travelers  are  found  in  the  particle  array  of 
the  destination  cell,  the  incoming  particle  is  simply  added  to  the  end. 

After  careful  optimization,  this  scheme  wa.s  found  to  run  nearly  as  fast  as  the 
communication  link  case,  but  with  a  20%  reduction  in  memory  requirements  for  a 
half-million  particle  calculation. 

A  code  listing  of  the  communication  function  is  included  in  Appendix  B. 

2.7  Boundary  Conditions 

As  for  any  solution  method,  proper  specification  of  the  boundary  conditions  is  critical 
to  a  successful  DSMC  simulation.  The  means  of  specifying  these  conditions  in  particle 
methods,  however,  differs  considerably  from  continuum  CFD.  In  continuum  CFD, 
macroscopic  variables,  such  as  temperature  and  velocity,  are  imposed  at  given  points 
based  on  their  intended  physical  state.  For  particle  methods,  these  conditions  must 
be  translated  into  rules  for  treating  individual  particles  near  these  points. 

2.7.1  Solid  Walls 

As  noted  for  intermolecular  collisions,  many  details  of  gas-surface  interaction  are  still 
unknown.  Again,  models  which  reproduce  the  important  features  of  the  physical 
situation  have  been  developed.  The  most  common  of  these  divides  particle  reflection 
from  solid  surfaces  into  two  classes:  specular  and  diffuse.  A  given  interaction  may  be 
described  completely  by  one  of  these  classes  or  by  some  combination  of  the  two.  This 
description  is  quantified  by  the  tangential  momentum  accommodation  coefficient, 
F,  which  varies  from  0,  for  no  accommodation  (specular  reflection),  to  1,  for  full 
accommodation  (diffuse  reflection).  In  the  code,  this  is  accomplished  by  treating  F 
as  the  probability  a  given  particle  reflection  will  be  treated  diffusely. 
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Specular  Reflection 

In  a  specular  reflection,  the  normal  component  of  the  impinging  particle’s  velocity 
vector  is  simply  reversed  and  the  tangential  component  is  left  unchanged.  No  modifi¬ 
cation  is  made  to  the  energy  of  the  particle.  This  is  intended  to  model  an  interaction 
with  a  perfectly  smooth  (ie.  frictionless)  surface.  It  may  also  be  used  to  simulate  a 
symmetry  plane. 

Diffuse  Reflection 

In  a  diffuse  reflection,  the  impinging  particle  is  emitted  from  the  surface  without  re¬ 
gard  to  its  incoming  state.  The  outgoing  velocity  is  randomly  assigned  according  to 
a  half-range  Maxwellian  distribution  at  the  wall  temperature.  This  is  known  as  full 
thermal  and  momentum  accommodation  and  may  be  viewed  in  terms  of  a  particle 
that  is  absorbed,  then  re-emitted  at  equilibrium  with  the  surface.  This  is  intended 
to  model  an  interaction  with  a  completely  rough  surface,  which  is  considered  a  valid 
description  of  most  engineering  materials.  MEMS  devices,  however,  often  contain 
surfaces  which  are  cut  along  the  crystal  planes  of  Silicon  wafers.  Tangential  momen¬ 
tum  accommodation  coefficients  considerably  less  than  one  are  therefore  possible  on 
these  extraordinarily  smooth  surfaces. 

2.7.2  Inflow/Outflow  Faces 

Inflow/outflow  (I/O)  faces  are  considerably  more  difficult  to  treat  than  solid  bound¬ 
aries,  particularly  for  low-speed  cases,  such  as  those  presented  in  Sections  3.1  and 
3.2.  Ironically,  this  task  appears  to  be  straightforward  and  well-defined;  simply  in¬ 
troduce  and  remove  particles  to  obtain  the  desired  flow  state.  Upon  closer  inspection, 
however,  it  is  found  to  be  another  sensitive  compromise  between  physical  detail  and 
computational  efficiency.  In  this  case,  an  improper  formulation  leads  to  a  strong 
flow  adjustment  (very  similar  to  a  shock)  near  the  boundary.  This  effect  is  shown 
in  Figure  2-3  for  a  channel  with  a  specified  pressure  ratio  of  4,  which  relaxed  to  ap¬ 
proximately  3.65.  Unfortunately,  the  strength  of  this  relaxation  is  difficult  to  predict, 
frustrating  efforts  to  model  a  specific  geometry  and  flow  condition. 

In  the  current  formulation,  I/O  faces  are  treated  in  two  stages  of  the  calculation; 
particle  movement  and  boundary  enforcement. 

During  the  movement  stage,  particles  are  simply  removed  from  the  calculation  if 
they  encounter  an  I/O  face.  The  data  structure  and  communication  schemes  outlined 
previously  make  this  a  straightforward  operation;  the  particle’s  position  in  its  cell  is 
marked  vacant,  as  if  it  has  moved  to  another  cell,  but  it  is  not  labelled  as  a  traveler, 
effectively  moving  it  ‘nowhere’. 

During  the  enforcement  stage,  particles  are  introduced  at  I/O  faces  to  maintain 
user-specified  boundary  conditions,  which  are  expressed  in  terms  of  macroscopic  vari¬ 
ables.  These  boundary  conditions,  combined  with  quantities  calculated  from  the 
current  flow  state,  determine  the  number  of  particles  to  introduce  and  their  velocity 
distribution,  as  detailed  below.  It  should  be  noted  that,  although  the  mechanism 
for  enforcing  boundary  conditions  differs  considerably  between  particle  methods  and 
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Figure  2-3:  Pressure  distribution  for  a  channel  with  poorly-formulated  10  treatment. 

continuum  CFD,  the  choice  of  which  macroscopic  variables  to  specify  externally  and 
which  to  calculate  from  the  domain  is  determined  in  both  techniques  by  the  ‘charac¬ 
teristic  lines’. 

Characteristics  lines,  or  simply  ‘characteristics’,  are  paths  in  space  and  time,  de¬ 
rived  from  the  Euler  equations,  along  which  certain  flow  variables  remain  constant. [20] 
They  are  therefore  said  to  “carry”  information  from  one  place  to  another.  Charac¬ 
teristics  are  used  when  formulating  boundary  conditions  to  prescribe  how  much  in¬ 
formation  is  communicated  to  the  boundary  from  inside  the  domain  and  how  much 
is  communicated  from  outside.  This  determines  which  variables  may  be  specified  and 
which  must  be  calculated  from  the  flow  itself. 

There  are  four  characteristic  lines;  one  carries  the  entropy,  one  carries  the  trans¬ 
verse  speed,  and  two  carry  the  Riemann  Invariants,  J+  and  J_,  which  are  given  by: 

T  .  2a 

J±=u± - -  (2.6) 

where  u  is  the  flow  speed,  a  is  the  speed  of  sound,  7  is  the  ratio  of  specific  heats,  and 
all  quantities  are  dimensional.  The  first  two  of  these  move  through  space  with  speed 
u  and  the  second  two  with  speed  u  +  a  and  u  -  a,  respectively.  Thus,  three  of  the 
characteristics  always  point  in  the  flow  direction  and  the  fourth  points  upstream  for 
subsonic  flow  {u  <  a)  and  downstream  for  supersonic  flow  (u  >  a).  A  subsonic  inlet 
therefore  takes  information  from  outside  the  domain  on  three  characteristics  and  from 
inside  the  domain  on  one.  A  supersonic  inlet,  however,  takes  all  of  its  information 
from  outside  the  domain.  Similarly,  a  subsonic  outlet  takes  three  characteristics 
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from  inside  the  domain  and  one  from  outside  while  a  supersonic  outlet  takes  all  its 
characteristics  from  inside  the  domain.  Over-constraining  the  boundary,  by  specifying 
too  many  variables  for  the  number  of  incoming  characteristics,  leads,  in  DSMC,  to 
strong  local  flow  adjustments  of  the  type  shown  in  Figure  2-3. 

Although  the  characteristics  constrain  the  number  of  state  variables  that  can 
be  specified  at  the  I/O  faces,  there  is  some  freedom  to  choose  their  identity.  [21]  One 
possible  arrangement  is  to  specify  the  streamwise  speed,  transverse  speed,  and  density, 
calculating  the  pressure  from  the  domain.  This  is  useful  for  modeling  aerodynamic 
bodies  at  a  given  flight  speed,  angle  of  attack,  and  altitude,  for  example.  A  more 
common  arrangement  is  obtained  by  substituting  temperature  for  density  in  the  above 
case  so  angle  of  attack  and  Mach  number  are  the  input  parameters.  Another  possible 
arrangement  was  used  for  the  calculations  presented  in  Chapter  3.  For  these  flows,  the 
pressure,  temperature,  and  transverse  speed  were  specified  and  the  streamwise  speed 
was  calculated  from  the  domain  at  inflow  faces.  Only  pressure  was  specified  at  the 
outflow  faces  of  the  subsonic  cases  and  nothing  was  specified  in  the  supersonic  cases. 
These  boundary  conditions  are  intended  to  model  fully-developed  flows  in  channels 
and  nozzles,  (ie.  flows  which  appear  to  represent  a  segment  of  a  device  which  is  far 
from  its  pressure  reservoirs) .  This  is  accomplished  by  enabling  the  streamwise  speed 
to  self-adjust  to  a  parabolic  profile  at  the  inlet  and  outlet  which  smoothly  blends  with 
the  velocity  profiles  in  the  remainder  of  the  channel. 

DSMC’s  statistical  nature  complicates  the  boundary  enforcement  process,  how¬ 
ever.  To  determine  a  cell’s  macroscopic  variables,  its  microscopic  state  must  be 
sampled.  Unfortunately,  there  may  be  as  few  as  20  particles  in  a  given  cell  at  a  given 
time,  so  a  single  sample  produces  unacceptable  statistical  scatter.  This  leaves  two 
options;  neighboring  cells  may  be  included  in  the  instantaneous  sample,  or  it  may  be 
replaced  with  a  time  average.  The  former  option  involves  the  selection  of  neighboring 
cells  with  states  ‘close-enough’  to  that  of  the  cell  in  question  to  yield  a  meaningful 
spatial  average.  The  latter  option  involves  choosing  a  time-averaging  method  that 
results  in  a  suflflciently  accurate  estimate  but  still  allows  the  flow  to  reach  its  steady 
state  with  reasonable  speed.  The  latter  option  was  implemented  in  the  current  code. 
The  cell  state  was  sampled  after  a  movement  step  was  completed,  incoming  particles 
were  introduced  at  boundary  enforcement,  and  collisions  were  performed.  A  weighted 
average  was  then  taken  between  this  result  and  a  running  value  collected  from  previ¬ 
ous  time  steps.  The  weighting  of  the  instantaneous  state  may  be  varied  to  balance  the 
accuracy  of  the  running  estimate  with  the  convergence  speed;  a  large  weight  makes 
the  estimate  sensitive  to  the  statistical  scatter  inherent  in  the  instantaneous  average, 
causing  error  in  the  estimate,  while  a  small  weight  causes  prior  information  to  decay 
slowly,  retarding  convergence  to  steady-state.  A  weight  of  1/20  was  chosen  for  the 
cases  presented  in  this  work. 

Once  the  necessary  macroscopic  variables  are  obtained  from  either  the  user  or  the 
domain,  the  boundary  enforcement  process  is  identical  for  inflow  and  outflow  faces: 
particles  are  introduced  in  sufficient  quantity  and  with  the  proper  velocity  distribution 
to  satisfy  local  constraints.  It  should  be  noted  that  a  significant  number  of  particles 
are  introduced  at  both  the  inflow  and  the  outflow  faces  in  low-speed  calculations, 
consistent  with  the  existence  of  the  backward-running  characteristic  discussed  above. 
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First,  the  number  of  particles  to  introduce  at  an  I/O  face  must  be  calculated.  To 
accomplish  this,  a  target  number  density  for  the  cell  containing  the  face  is  obtained 
from  the  specified/calculated  macroscopic  variables  and  the  ideal  gas  relation,  which 
is  simply 

p  =  nT  (2.7) 

under  the  non-dimensionalizations  of  Section  2.1.  A  target  particle  quantity  is  then 
calculated  through  multiplication  by  the  cell  volume.  The  actual  particle  quantity 
is  compared  to  this  target  to  determine  the  number  of  particles  to  introduce.  If  the 
actual  particle  quantity  is  greater  than  or  equal  to  the  target,  no  action  is  taken; 
particles  are  never  removed  from  a  cell. 

Velocity  components  perpendicular  to  the  I/O  face  are  then  assigned  to  the  in¬ 
coming  particles  according  to  a  Maxw^ellian  distribution  at  the  specified/calculated 
temperature.  In  non-dimensional  form,  the  velocity  distribution  for  a  thermal  velocity 
component  in  the  transverse  direction,  say  v',  is  given  by: 

=  (2.8) 

Values  are  selected  from  this  distribution  directly  using  a  method  presented  by  Bird 
in  Ref.  [1]:  select  two  independent  random  numbers,  Rf^  and  Rf^,  then  calculate  v' 
from: 


6  =  277  Rf^ 

r  ^  ^-rin(%)  (2.9) 

v'  =  rsm{0) 

At  inflow  faces,  the  specified  transverse  velocity  (which  is  zero  for  all  cases  in  this 
report)  is  added  to  this  value.  At  outflow  faces,  the  calculated  transverse  velocity 
is  added.  This  operation  is  identical  for  all  cases  because  one  of  the  characteristics 
discussed  above  carries  the  transverse  velocity  specifically  and  always  moves  in  the 
flow  direction. 

The  velocity  component  normal  to  the  I/O  face  is  then  assigned  according  to 
a  fluxal  distribution,  ie.  a  distribution  which  is  shifted  based  on  the  mean  normal 
velocity  through  the  face.  The  resulting  (non-normalized)  velocity  distribution  is 
given  by: 

A-  (2.10) 

where  Un  is  the  normal  component  of  the  macroscopic  (mean)  velocity  across  the  face 
and  a  positive  value  denotes  motion  into  the  cell  (for  both  u'  and  rt„).  This  distri¬ 
bution  is  sampled  via  the  acceptance-rejection  method,  where  a  randomly-selected 
value  for  u'  is  either  accepted,  with  a  probability  equal  to  fu'{u')/{fu')max,  or  rejected 
and  another  value  chosen,  repeating  until  an  acceptance  is  made.  Situations  where 
u'  <  0  are  non-physical,  because  the  particle  could  not  enter  the  cell,  and  are  therefore 
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excluded  from  consideration. 

A  low  aspect  ratio  channel  with  a  pressure  ratio  of  3  and  the  fully-developed  flow 
boundary  conditions  described  above  served  as  the  test  case  during  10  boundary 
development.  This  case  was  chosen  because  it  produces  strong  gradients  yet  remains 
subsonic  at  the  outflow,  presenting  a  significant  challenge  to  boundary  enforcement. 
The  channel  geometry  was  150x30  A  (referenced  to  the  inlet)  with  4000  uniform  cells 
and  approximately  180,000  particles  at  steady-state.  A  sample  run  is  presented  to 
demonstrate  the  effectiveness  of  the  current  boundary  formulation. 

Figure  2-4  contains  the  complete  (ie.  all  cells  were  plotted)  streamwise  velocity 
distribution.  End-effects  are  almost  imperceptable  at  both  the  inlet  and  the  outlet. 
The  “long  channel”  boundary  conditions  were  therefore  successful  and  the  flow  in 
the  entire  domain  can  be  considered  fully-developed.  The  corresponding  pressure 
distribution  is  shown  in  Figure  2-5,  also  with  good  results,  exhibiting  very  little  of 
the  flow  adjustment  present  in  Figure  2-3. 


Figure  2-4:  Complete  streamwise  velocity  distribution  for  I/O  test  channel. 

A  code  listing  of  the  boundary  enforcement  function  is  given  in  Appendix  C. 

2.8  Verification 

Before  treating  more  complicated  cases,  simple  verification  runs  were  made  to  test 
the  algorithm  formulation  and  coding.  Cases  with  either  an  analytical  solution  or 
published  results  were  chosen  for  this  task  to  facilitate  comparison  with  the  current 
code.  An  equilibrium  gas  in  a  box  and  a  1-D  shock  wave  are  examined  below. 
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Figure  2-5:  Complete  pressure  distribution  for  I/O  test  channel. 

2.8.1  Equilibrium  Gas 

All  flows  examined  in  this  work  have  Knudsen  numbers  which  place  them  comfortably 
in  the  collisional  regime.  The  proper  treatment  of  inter-particle  encounters  is  therefore 
crucial  to  an  accurate  simulation.  To  test  the  mechanisms  governing  the  collision 
frequency,  the  equilibrium  collision  rate  of  several  different  gases  was  calculated  by 
simulating  a  resting  fluid  in  a  closed  domain.  The  results  were  then  compared  to  the 
theoretical  prediction. 

The  theoretical  equilibrium  particulate  collision  rate,  Vp,  represents  the  average 
number  of  collisions  suffered  per  particle  per  unit  time.  It  is  given  by  the  mean 
molecular  speed,  c',  divided  by  the  mean  free  path,  A.  From  Section  2.1,  the  non- 
dimensionalization  factors  for  length  and  speed  are  A  and  c^,  respectively,  both  at 
a  reference  condition.  Noting  that  c'  =  2/y/K  the  non-dimensional  equilibrium 
collision  rate  is  found  to  be: 


j/p  =  4=  =  1.13  (2.11) 

VTT 

To  verify  that  the  DSMC  code  correctly  produces  this  rate,  a  10x10  domain, 
composed  of  100  cells  and  specularly-reflecting  walls,  was  constructed.  Five-thousand 
particles  were  distributed  in  this  domain  with  velocities  selected  from  a  Maxwellian 
distribution  at  the  reference  temperature.  The  flow  was  then  run  for  1000  time  steps, 
to  allow  any  initial  condition  effects  to  die  out,  then  sampled  every  four  time  steps 
until  1000  samples  were  obtained.  This  was  performed  for  four  gases  with  significantly 
different  values  of  uj  and  the  results  are  presented  in  Table  2.2. 

First,  it  may  be  noted  that  the  collision  rate  is  independent  of  molecular  species. 
This  is  supported  by  the  theoretical  prediction,  as  no  gas  constants  appear  in  Equa¬ 
tion  2.11.  In  addition,  the  computed  rate  is  within  2%  of  the  theoretical  value.  This 
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Gas 

LU 

lyp 

Neon 

0.16 

1.15 

Nitrogen 

0.24 

1.15 

Argon 

0.31 

1.15 

Carbon  Dioxide 

0.43 

1.15 

Table  2.2:  Computed  Equilibrium  Collision  Rates  for  Various  Gases 


error  was  found  to  be  a  function  of  the  number  of  particles  per  cell, 
dance  with  the  statements  made  in  Section  1.4.1.  This  dependence  is 
in  Table  2.3. 


Nc,  in  accor- 
demonstrated 


Nc 

Up 

%Error 

100 

1.14 

1.0 

50 

1.15 

1.9 

20 

1.17 

3.7 

10 

1.23 

9.0 

5 

1.32 

17.0 

Table  2.3:  Computed  Equilibrium  Collision  Rates  for  Various  Cell  Particle  Numbers 


2.8.2  1-D  Shock  Wave 

To  assess  the  accuracy  of  the  non-equilibrium  collision  rate  and  the  post-collision 
scattering  law,  a  1-D  shock  wave  was  computed.  This  case  was  selected  because 
analytical  expressions  exist  for  the  state  variable  jumps  across  the  shock  and  published 
DSMC  results  for  the  shock  profile  are  available. 

This  calculation  was  performed  on  a  100x40  grid  with  4000  uniform  rectangular 
cells.  The  simplicity  of  this  geometry  allowed  a  particle’s  current  cell  to  be  determined 
by  mathematical  means,  based  on  its  position,  the  grid  dimensions,  and  the  number  of 
cells.  The  trajectory-tracing  movement  scheme  outlined  in  Section  2.4  was  therefore 
not  required.  This  enabled  the  shock  to  be  created  by  moving  the  left  side  of  the 
domain,  compressing  the  entire  grid  with  each  time  step  to  maintain  its  regularity. 
The  grid  was  initialized  with  a  resting  fluid,  and  the  wall  was  set  in  motion  at  t=0. 
Ensemble  averaging  was  performed  for  a  total  of  50  ensembles. 

The  resulting  jump  in  state  variables  across  the  shock  may  be  compared  to  the 
analytical  predictions  derived  from  the  1-D  flow  equations.  [20]  This  comparison  is 
presented  in  Table  2.8.2  for  a  Mach  8  shock  in  Argon  {u^aii  =  5.39,  u  =  0.31).  An 
excellent  agreement  is  obtained  in  all  cases,  with  a  maximum  error  of  less  than  2%. 

The  spatial  profile  of  the  shock  may  also  be  compared  to  other  DSMC  results, 
such  as  those  presented  by  Baganoff  and  McDonald.  [11]  This  comparison  is  shown  in 
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Variable 

Analytical 

Computed 

%Error 

Pressure 

79.75 

79.54 

-0.26 

Temperature 

20.87 

20.58 

-1.39 

Density 

3.82 

3.86 

-hi. 05 

Wave  Speed 

8.00 

8.08 

-hi. 00 

Table  2.4:  Analytical  and  Computed  State  Variable  Jumps  Across  a  Mach  8  Shock 


Figure  2-6  for  a  Mach  3  shock  in  a  Maxwellian  gas,  which  is  equivalent  to  a  VHS  gas 
with  a;  =  0.5.  The  agreement  is,  again,  excellent,  providing  evidence  that  the  code  is 
indeed  properly  simulating  the  fluid  behavior. 


Figure  2-6:  Comparison  of  Mach  3  shock  profiles  from  the  current  code  and  published 
DSMC  results. 
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Chapter  3 

Non- Continuum  Results 


As  discussed  in  Chapter  1,  DSMC  is  primarily  used  for  flows  in  non-continuum 
regimes,  for  which  Navier-Stokes  based  methods  break  down.  Formerly,  most  of 
these  flows  involved  high-altitude  flight.  Now,  MEMS  devices  provide  a  rich  array  of 
interesting  flows  in  these  regimes  that  have  practical  applications  in  a  wide  variety  of 
areas.  A  sampling  of  these  cases  are  examined  in  this  chapter.  This  effort  is  intended 
to  investigate  DSMC’s  ability  to  accurately  and  efficiently  model  micro-flows  as  well 
as  illuminate  some  of  their  unique  features  which  are  important  to  MEMS  designers. 


3.1  Slip  Flow  Regime  Micro-Channel 

The  first  case  explored  was  a  steady  flow  through  a  micro-channel  with  an  outlet 
Knudsen  number  of  0.05,  placing  it  in  the  slip  flow  regime.  This  geometry  is  similar 
to  those  investigated  experimentally  by  Harley  et  a/.  [22]  and  Arkilic  et  a^.[23]  and 
numerically  (with  spectral  element  methods)  by  Beskok  and  Karniadakis.[24]  This  is, 
historically,  an  important  canonical  case  for  determining  the  effect  of  rarefaction  on 
the  transport  terms  in  the  Navier-Stokes  equations.  It  is  also  useful  as  a  representation 
of  the  flow  along  certain  features  common  in  MEMS  devices,  such  as  the  space  under 
the  floating  plate  of  a  shear  stress  sensor  or  accelerometer,  which  is  typically  only  one 
micron  high  but  hundreds  of  microns  in  breadth  and  depth.  In  addition,  it  serves 
as  an  interesting  calibration  case  to  assess  the  accuracy  of  the  numerical  algorithm 
because  analytical  solutions  have  been  developed  for  this  geometry. 

In  one  such  effort,  Arkilic  et  al.  show  that  the  Navier-Stokes  equations  may  be 
solved  analytically  for  a  long,  high  aspect-ratio,  isothermal  channel  in  the  slip  flow 
regime  if  the  boundary  conditions  are  modified  to  include  a  ATn-dependent  streamwise 
velocity  (slip)  at  the  wall,  given  by: 


fJ'watl 


2-F 

F 


Kn 


(3.1) 


where  u  is  the  streamwise  velocity,  Kn  is  the  local  Knudsen  number,  y  is  the  transverse 
coordinate,  which  has  its  zero  at  the  channel  centerline,  and  F  is  the  tangential 
momentum  accommodation  coefficient,  discussed  in  Section  2.7.1. 
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Through  this  analysis,  an  expression  may  be  obtained  for  the  pressure  distribution 
in  a  microchannel  with  diffusely-reflecting  walls  as  a  function  of  streamwise  channel 
location  and  overall  pressure  ratio: 

V{x)  =  -6Kno  +  |  [{Vf  -  1)  +  l2Kno{V^  -  1)]  (3.2) 

where  V{x)  and  are  the  local  and  inlet  pressures,  respectively,  normalized  by  the 
outlet  value,  Krio  is  the  outlet  Knudsen  number,  x  is  the  streamwise  coordinate, 
and  L  is  the  channel  length.  The  distribution  predicted  by  Equation  3.2  may  be 
compared  to  a  DSMC  result  as  a  test  of  both  theory  and  code.  Such  a  comparison  is 
presented  in  Figure  3-1  for  a  600  x  20A  (referenced  to  the  outlet)  channel  run  with 
Nitrogen  {u  =  0.24)  at  a  pressure  ratio  of  2.47  with  the  infinite  channel  boundary 
conditions  discussed  in  Section  2.7.2.  A  good  agreement  is  obtained  (max  error  = 
1.5%),  including  the  nonlinear  pressure  distribution  that  occurs  due  to  the  large 
pressure  drop  down  the  length  of  the  channel. 


Figure  3-1:  Comparison  of  computed  and  analytical  pressure  distributions  for  a  micro- 
channel  in  the  slip  flow  regime. 

A  theoretical  expression  for  the  streamwise  velocity  distribution  was  also  devel¬ 
oped  by  Arkilic  et  ai: 


/  2 

2/i  dx  \ 


(3.3) 


where  /r  is  the  coefficient  of  viscosity,  p  is  the  pressure,  and  H  is  the  channel  height. 
This  equation  is  plotted  in  Figure  3-2  for  the  geometry  and  conditions  used  above. 
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Several  features  unique  to  a  flow  of  this  type  are  visible.  First,  the  fluid  accelerates 
as  it  moves  down  the  channel,  unlike  in  the  familiar  Poiseuille  result.  This  is  a 
consequence  of  the  density  drop  caused  by  the  decreasing  pressure  in  the  streamwise 
direction  (the  flow  is  effectively  isothermal).  The  mean  streamwise  velocity  must 
therefore  increase  to  maintain  a  constant  mass  flow.  Second,  the  velocity  at  the  walls 
is  nonzero  and  increases  with  increasing  x-coordinate.  This  is  the  aforementioned 
‘slip  flow’,  which,  by  Equation  3.1,  is  essentially  zero  for  continuum  flows  due  to  their 
very  small  Knudsen  number.  The  increase  in  slip  velocity  down  the  channel  is  a  result 
of  growth  in  both  Kn  (from  the  decreasing  pressure)  and  velocity  gradient  at  the  wall 
(from  the  accelerating  flow). 


Figure  3-2:  Theoretical  streamwise  velocity  distribution  for  a  micro-channel  in  the 
slip-flow  regime. 

The  DSMC  result  for  this  configuration  is  presented  in  Figure  3-3.  Comparing 
this  to  the  previous  figure,  it  may  be  concluded  that  the  DSMC  calculation  qualita¬ 
tively  reproduces  the  mean  flow  acceleration  and  the  increasing  slip  flow  predicted 
by  the  theory.  In  addition,  good  quantitative  agreement  is  obtained  in  the  velocity 
distributions. 

A  further  comparison  with  the  theoretical  analysis  of  Arkilic  et  al.  may  be  found 
by  normalizing  the  velocity  distribution  of  Equation  3.3  by  the  average  velocity  at  a 
given  x-location,  obtained  by  integrating  u  from  the  lower  wall  to  the  upper  wall  and 
dividing  by  H\ 
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Figure  3-3:  Computed  streamwise  velocity  distribution  for  the  micro-channel  of  Fig¬ 
ure  3-2. 


SI-« 


Rearranging  the  resulting  expression  yields: 


The  left  side  of  this  equation,  which  will  be  referred  to  as  the  ‘similarity  speed’, 
Us,  is  a  function  of  x  and  y,  while  the  right  is  a  function  of  y  only.  Consequently,  if 
the  slip-flow  analysis  holds,  calculating  the  similarity  speed  using  the  local  Kn{x)  and 
u{x,  y)  will  yield  identical  parabolas  at  all  rc-stations  down  the  length  of  the  channel. 

This  assertion  was  tested  by  computing  a  similarity  speed  distribution  from  the 
DSMC  output.  The  maximum  and  minimum  of  the  result  at  each  x-location  is  shown 
in  Figure  3-4.  It  should  be  noted  that  the  upper  theoretical  line  is  not  placed  exactly 
at  0.25  because  an  even  number  of  cells  was  used  in  the  DSMC  run,  so  there  is  no 
data  point  in  the  center  of  the  channel.  Also,  the  lower  line  is  not  at  zero  because 
the  macroscopic  quantities  for  a  cell  are  assumed  to  be  associated  with  its  centroid, 
so  there  are  no  data  points  on  the  walls  themselves. 

It  may  be  concluded  from  this  figure  that  the  similarity  assertion  indeed  holds  for 
the  slip  flow  channel.  As  predicted  by  the  analysis,  the  down-channel  variation  of  the 
streamwise  velocity  profile  seen  in  Figure  3-3  has  given  way  to  a  constant  similarity 
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Figure  3-4:  Comparison  of  computed  and  theoretical  maximum  and  minimum  simi¬ 
larity  speeds  for  a  micro-channel  in  the  slip  flow  regime. 

speed  profile.  In  addition,  the  maximum  and  minimum  similarity  speeds  compare 
well  with  the  predicted  values. 

Overall,  excellent  agreement  was  obtained  between  the  analytical  solution  of  Ark- 
ilic  et  al.  and  the  DSMC  results.  This  supports  the  accuracy  of  both  techniques. 
For  the  DSMC  code,  however,  it  is  just  the  beginning;  many  more  interesting  flows, 
for  which  there  are  no  reliable  analytical  solutions,  may  be  easily  treated  with  this 
method.  The  remaining  cases  presented  in  this  chapter  are  intended  to  demonstrate 
this  capability. 

3.2  Transition  Regime  Micro-Channel 

One  of  the  great  strengths  of  DSMC  is  its  validity  for  dilute  gases  in  all  Knudsen 
number  regimes.  One  of  the  most  interesting  of  these  is  the  transition  regime,  de¬ 
fined  in  Section  1.2  as  0.1  <  Kn  <  3.  Here  the  mean  free  path  is  comparable  to 
the  characteristic  dimension  of  the  flow.  This  makes  analytical  solution  very  difficult 
because  the  approximation  of  transport  terms  based  on  macroscopic  quantities  be¬ 
comes  inaccurate,  precluding  the  use  of  the  Navier-Stokes  equations  (even  with  the 
slip-flow  boundary  condition  of  Equation  3.1).  Collisions  are  still  important,  however, 
so  the  collisionless  Boltzmann  equation  is  not  yet  an  option.  This  leaves  only  the  full 
Boltzmann  equation;  a  very  difficult  expression  to  solve,  either  analytically  or  with 
numerical  techniques. 


33 


DSMC  is  therefore  a  very  attractive  tool  for  investigating  the  transition  regime.  It 
IS  also  one  of  the  few  tractable  techniques  which  is  uniformly  valid  in  mixed  Kn-regime 
flows.  These  are  important  features  because,  due  to  the  aforementioned  difficulties, 
relatively  little  is  known  about  these  cases.  Such  knowledge  is  critical,  however, 
because  many  MEMS  devices  contain  flows  of  this  nature. 

The  micio-channel  was  again  used  as  a  sample  case,  this  time  to  observe  the  failure 
of  the  slip-flow  analysis  as  Kn  enters  the  transition  regime.  The  channel  treated  here 
is  136.4  X  2.3  A,  with  a  pressure  ratio  of  4.2  and  an  outlet  Knudsen  number  of  0.44. 
The  working  fluid  is,  again.  Nitrogen. 


Figure  3  5:  Comparison  of  computed,  slip-flow,  and  continuum  pressure  distributions 
for  a  micro-channel  in  the  transition  regime. 

Proceeding  as  in  the  slip-flow  case,  the  computed  pressure  distribution  is  compared 
to  the  slip-flow  prediction  in  Figure  3-5.  The  continuum  curve  {Kn  =  0.0)  is  also 
shown  for  reference.  As  expected,  the  excellent  agreement  between  DSMC  and  theory 
obtained  for  the  slip-flow  case  (Figure  3-1)  is  no  longer  present.  The  error  between  the 
curves  has  grown  from  less  than  2%  to  more  than  4%.  The  form  of  this  disagreement 
is  also  significant:  the  computed  curve  is  more  linear  than  its  analytical  counterpart. 
A  trend  of  increasing  pressure  curve  linearity  with  increasing  rarefaction  is  therefore 
established  by  the  relative  shape  of  the  continuum,  slip-flow,  and  transition  curves. 

The  analytical  prediction  for  the  streamwise  velocity  distribution  in  this  channel 
IS  presented  in  Figure  3-6.  Note  that  the  theory  predicts  a  flatter  profile  than  for 
the  previous  case.  This  may  be  attributed  to  the  last  term  in  Equation  3.3,  which 
is  constant  across  the  channel  and  proportional  to  Kn.  The  Knudsen  number  is  an 
order  of  magnitude  larger  in  this  case,  so  this  term  has  a  much  stronger  influence  on 
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the  shape  of  the  distribution. 


Figure  3-6:  Theoretical  streamwise  velocity  distribution  for  a  transition-regime  micro 
channel  calculated  with  the  slip-flow  analysis. 

Upon  plotting  the  computed  distribution  for  comparison  (Figure  3-7),  it  becomes 
evident  that  the  assumptions  supporting  Equation  3.3  are  beginning  to  fail.  Both 
the  slip  flow  and  maximum  speeds  at  a  given  x-location  are  higher  than  predicted  by 
as  much  as  40%.  This  allows  the  channel  to  support  a  much  larger  mass  flow  than 
would  be  predicted  by  the  slip  flow  theory,  which,  in  turn,  predicts  a  larger  mass  flow 
than  the  traditional  Navier-Stokes  analysis. 

A  final  exposition  of  transition  behavior  may  be  made  via  the  similarity  analysis  of 
the  previous  section.  As  discussed  in  Section  1.2,  the  transition  regime  is  commonly 
considered  to  begin  at  Kn  =  0.1.  This  assertion  may  be  tested  by  noting  that  Kn 
increases  with  downstream  position.  The  similarity  profiles,  which  depend  on  the  slip- 
flow  solution,  may  then  be  computed  for  each  position  and  compared  to  the  analytical 
prediction.  Because  each  position  has  a  Knudsen  number  associated  with  it,  the  value 
for  which  the  slip  flow  analysis  fails  may  be  determined  by  finding  the  point  where 
the  experimental  and  analytical  curves  begin  to  diverge.  Toward  this  end.  Figure  3-8 
contains  the  computed  maximum  and  minimum  similarity  speeds,  plotted  with  the 
analytical  prediction  in  a  fashion  identical  to  that  of  Figure  3-4.  The  Kn  distribution 
was  then  overlaid  to  facilitate  determining  its  value  when  the  slip-flow  analysis  fails. 

It  may  be  concluded  from  this  figure  that  the  slip  flow  analysis  begins  to  fail  at 
approximately  Kn  =  0.15.  This  supports  the  oft-used  boundary  for  the  transition 
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Figure  3-7:  Computed  streamwise  velocity  distribution  for  the  channel  of  Figure  3-6 

region,  Kn  =  0.1.  This  limit  may  be  understood  if  the  slip  boundary  condition, 
Equation  3.1,  is  viewed  as  an  expansion  of  the  wall  velocity  in  powers  of  Kn.  The 
no-slip  condition  is  then  the  zeroth-order  solution,  and  Equation  3.1  is  the  first-order 
solution.  It  is  therefore  logical  that  the  neglected  higher-order  terms  would  begin  to 
significantly  affect  the  result  when  Kn  exceeds  0.1.  A  second-order  accurate  boundary 
condition  in  terms  of  the  continuum  variables  is  presented  in  Ref.  [24],  though  it 
should  be  remembered  that  the  Navier-Stokes  equations  themselves  are  only  strictly 
valid  to  first  order  in  Kn. 

3.3  Supersonic  Micro-Nozzles 

The  final  cases  presented  in  this  chapter  are  supersonic  micro-nozzles.  These  may 
be  viewed  as  channels  whose  upper  and  lower  walls  form  a  parabolic  contraction  / 
expansion.  Two  such  nozzles  are  discussed,  one  with  an  area  ratio  of  3.5  and  sonic 
flow  at  the  throat,  and  one  with  an  area  ratio  of  2.0  and  subsonic  flow  at  the  throat. 
The  grid  for  the  latter  case  is  shown  in  Figure  3-9. 

Analytically,  a  sufficiently  long  nozzle  can  be  considered  quasi- ID  and  the  solu¬ 
tion  given  in  Section  3.1  may  be  used,  with  appropriate  modifications  for  the  slowly 
varying  channel  height.  The  nozzles  presented  in  this  section,  however,  have  a  total 
length  of  only  six  times  their  throat  height,  so  the  quasi-lD  assumption  is  not  valid. 
In  addition,  the  significant  expansion  may  cause  the  Ah-regime  to  change  at  one  or 
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Figure  3-8:  Computed  and  analytical  similarity  speeds  with  Knudsen  number  overlay. 

more  streamwise  positions.  These  factors  make  analytical  and  continuum-based  nu¬ 
merical  treatment  of  these  geometries  difficult.  Nonetheless,  nozzles  such  as  these 
may  play  important  roles  in  devices  such  as  micro-rocket  thrusters  and  micro-gas 
turbine  generators,  for  example.  Investigating  their  behavior  is  therefore  a  valuable 
task  for  which  DSMC  is  well-suited. 

Both  cases  were  run  with  an  inlet  at  10  times  the  reference  pressure  and  ex¬ 
hausted  to  ‘vacuum’.  The  latter  condition  was  implemented  by  simply  removing 
from  the  simulation  any  particle  which  crossed  the  outlet  boundary  and  refraining 
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Figure  3-9:  Computational  grid  for  a  micro-nozzle  with  area  ratio  2.0. 
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from  introducing  new  particles  at  those  faces.  The  walls  were  isothermal  at  the  ref¬ 
erence  temperature  with  full  thermal  and  momentum  accommodation.  The  working 
fluid  was  Helium  (a;  =  0.20),  which  was  supplied  at  the  reference  temperature. 

3.3.1  Sonic  Throat 

The  first  nozzle  presented  has  an  area  ratio  of  3.5  and  4200  cells.  With  the  boundary 
conditions  given  above,  the  resulting  pressure  ratio  was  approximately  24  and  the 
outlet  Knudsen  number,  based  on  the  passage  height,  was  0.03. 

The  Mach  number  distribution  for  this  case  is  shown  in  Figure  3-10.  A  number 
of  interesting  features  are  visible.  First,  as  normally  expected,  the  flow  is  sonic  at 
the  throat.  Second,  a  Mach  number  of  2.4  is  reached.  This  is  considerable  when  it 
is  noted  that  the  nozzle  is  only  approximately  600  inlet  mean  free  paths  in  length. 
Finally,  the  slip  flow  speed  is  substantial,  exceeding  Mach  0.5  near  the  outlet. 
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Figure  3-10:  Mach  number  distribution  for  a  micro-nozzle  with  a  sonic  throat. 

The  temperature  distribution  for  this  micro- nozzle  is  shown  in  Figure  3-11.  It 
is  clear  that,  unlike  the  channel  case,  this  flow  cannot  be  considered  isothermal.  A 
strong  temperature  gradient  exists  in  both  the  streamwise  and  transverse  directions, 
creating  another  obstacle  to  analytical  treatment.  This  is  also  a  very  notable  feature 
when  the  diminutive  dimensions  of  the  nozzle  are  considered.  Though  the  walls  are 
only  about  30  local  mean  free  paths  apart  at  the  exit  and  are  isothermal  with  full 
energy  accommodation,  the  fluid  is  still  able  to  realize  a  substantial  reduction  in 
temperature.  The  effect  of  rarefaction  has  therefore  been  to  significantly  reduce  the 
thermal  communication  between  the  wall  and  the  fluid.  This  assertion  is  supported 
by  noting  the  large  thermal  slip  at  the  wall. 
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Figure  3-11:  Temperature  distribution  for  a  micro-nozzle  with  a  sonic  throat. 

3.3.2  Subsonic  Throat 

The  second  nozzle  presented  is  identical  to  the  first,  except  its  area  ratio  is  reduced 
to  2.0  and  its  grid  to  2400  cells.  With  the  boundary  conditions  given  above,  the 
resulting  pressure  ratio  was  10.2  and  the  outlet  Knudsen  number  was  0.03. 

The  Mach  number  distribution  for  this  case  is  shown  in  Figure  3-12.  This  dis¬ 
tribution  was  plotted  in  the  same  manner  as  the  previous  case,  only  the  viewpoint 
was  shifted  to  look  along  the  y-axis.  An  interesting  feature  is  now  visible:  the  outlet 
flow  is  supersonic,  but  the  sonic  point  is  downstream  of  the  throat.  This  result  is  at 
odds  with  the  inviscid,  quasi- ID  conclusion  that  sonic  flow  may  only  be  attained  at 
a  point  of  minimum  area.  [20]  The  highly  viscous  nature  of  this  flow  (due  to  the  close 
proximity  of  the  walls)  invalidates  this  prediction,  however,  causing  the  gas  to  con¬ 
tinue  to  accelerate  downstream  of  the  throat  despite  the  fact  that  it  is  still  subsonic 
and  the  duct  is  diverging. 

Because  a  significant  portion  of  this  nozzle  is  subject  to  the  competing  effects 
of  deceleration  due  to  geometry  and  acceleration  due  to  viscosity,  its  outlet  Mach 
number  is  considerably  («  30%)  smaller  than  the  previous  case  for  similar  boundary 
conditions.  This  highlights  the  importance  of  proper  design  in  such  a  device.  As  one 
of  the  few  analysis  tools  valid  for  these  flows,  DSMC  has  considerable  value  to  such 
an  effort. 
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Chapter  4 
Scaling  Issues^ 


The  previous  chapter  demonstrated  that  a  conventional  DSMC  code  can  accurately 
treat  many  geometries  of  interest  to  MEMS  designers.  Unfortunately,  the  demanding 
cell  size  and  time  step  constraints  outlined  in  Chapter  1  make  DSMC  modeling  of 
even  certain  micro-devices  very  difficult.  For  example,  the  mean  free  path  and  most 
probable  molecular  velocity  of  air  at  standard,  sea  level  conditions  are  approximately 
60  nm  and  414  m/s,  respectively.  Sizing  the  cell  lengths  at  one  A  and  the  time  steps 
at  as  called  for  in  Chapter  1,  a  2-D  simulation  of  the  smallest  micro-channel  of 
Arkilic  et  al.  {1.33x5000  /rm)  [23]  would  require  over  a  million  cells  and,  consequently, 
at  least  ten  million  particles.  In  addition,  the  time  steps  would  be  only  36  ps,  causing 
unsteady  cases,  even  those  with  microsecond  time  scales,  to  be  very  expensive. 

In  response  to  this  situation,  the  consequences  of  relaxing  the  cell  size  and  time 
step  requirements  commonly  placed  on  DSMC  are  explored  in  this  chapter.  The  goal 
of  this  effort  is  to  allow  the  cells  in  these  calculations  to  be  sized  by  the  flow  gradients, 
as  in  continuum  CFD,  rather  than  by  the  molecular  scales.  This  is  a  complex  issue, 
however,  because  this  sizing  violates  some  of  DSMC’s  underlying  assumptions  as 
the  scale  length  of  the  gradients  increases,  reducing  the  Knudsen  number.  Certain 
features  of  the  physical  situation  are  therefore  altered  or  lost.  The  key  to  a  successful 
simulation,  and  the  goal  of  this  chapter,  is  to  identify  these  features  and  assess  their 
importance  to  the  flow. 


4.1  Scaling  Rules 


As  discussed  in  Section  1.4.3,  the  cell  size  in  a  DSMC  calculation  should  be  comparable 
to  the  mean-free-path  of  the  gas  so  collision  partners  may  be  selected  without  regard 
to  their  position  in  the  cell.  In  this  chapter,  it  will  be  convenient  to  express  this 
constraint  in  terms  of  a  ‘cell  Knudsen  number’,  Kuc' 


Kric  =  ~ —  >  1 
- 


(4.1) 


^The  author  is  indebted  to  David  Gonzales,  a  co-author  of  the  paper  (AIAA-95-2088)  on  which 
this  chapter  was  based,  for  allowing  the  material  to  be  presented  here. 
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where  Axrf  is  a  typical  cell  dimension  (again,  the  subscript  “d”  refers  to  dimensional 
quantities). 

A  second  requirement  (discussed  in  Section  1.4.2)  is  that  the  time  step,  be 
small  compared  to  a  characteristic  time  so  particle  movement  and  collisions  may  be 
decoupled; 


Atd  < 


(4.2) 


where  c'^  is  the  most  probable  molecular  speed. 

Borrowing  from  traditional  computational  mechanics,  this  constraint  may  be  con¬ 
veniently  expressed  as  a  ‘CFL’,  or  Courant-Friedrich-Lewy  condition: 


Axd 


<  1. 


(4.3) 


Physically,  satisfying  this  condition  requires  a  particle  to  reside  in  the  same  cell  for 
a  few  time-steps  to  provide  it  with  ample  opportunity  to  interact  with  other  particles. 
This  ensures  that  its  information  can  be  distributed  properly  through  the  computa¬ 
tional  domain.  It  is  important  to  realize  that,  unlike  its  CFD  counterpart,  this  CFL 
condition  is  not  a  stability  requirement,  but  a  validity  requirement.  Violation  of  the 
condition  will  still  yield  results,  but  they  may  be  inaccurate.  A  DSMC  computation 
must  therefore  be  constructed  with  very  careful  attention  to  its  scaling  constraints 
because  the  algorithm  itself  will  give  no  indication  that  its  results  are  completely 
non-physical. 

Under  the  NTC  collision  method  outlined  in  Section  1.4.3,  the  number  of  collision 
pairs,  Np,  to  be  considered  in  a  given  cell  is  computed  through  the  equation: 


Np  oc  Ai-7  (4.4) 

nVc 

where  Vc  is  the  normalized  cell  volume,  Nc  is  the  number  of  computational  particles  in 
the  cell,  and  h  is  the  reference  simulation  number  density,  which  is  the  total  number 
of  computational  particles  divided  by  the  non-dimensional  domain  volume. 

Noting  that  the  normalized  cell  volume  scales  like  Knc~^  and  the  simulation  num¬ 
ber  density  scales  like  NcKuc  ■,  it  may  be  concluded  that: 

^  oc  At.  (4.5) 

Thus,  for  a  fixed  size  computation  {Nc  held  constant),  Np  depends  only  on  the  time 
step  chosen  to  advance  the  solution. 

The  ratio  of  potential  collision  pairs  to  the  number  of  particles  in  a  cell  will  be 
referred  to  as  the  ‘over-collision  ratio’.  A  value  of  one  implies  that,  on  average,  each 
particle  will  considered  for  a  collision  during  every  time  step.  In  a  well-resolved 
DSMC  computation,  where  At  ss  0.2,  roughly  20%  of  a  given  cell’s  particles  will  be 
considered  for  collisions. 

The  preceding  expressions  reveal  the  essential  scaling  issues  involved  in  applying 
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DSMC  methods  to  \ow-Kn  flows:  as  the  physical  size  of  the  problem  increases,  it 
becomes  computationally  impractical  to  maintain  a  cell  size  on  the  order  of  the  mean- 
free-path.  It  must  therefore  increase  or,  using  the  terminology  introduced  above,  Kric 
must  decrease,  violating  the  constraint  of  Equation  4.1.  In  addition,  as  the  cell  size 
increases,  the  time  step  should  be  reconsidered.  Two  options  exist:  scale  Af  so  the 
CEL  number  remains  constant,  or  hold  it  at  some  small  value  and  let  the  CEL  number 
decrease. 

The  first  option  maintains  proper  information  propagation  across  the  domain, 
preserving  computational  efficiency.  Unfortunately,  increasing  At  also  causes  rapid 
growth  in  the  number  of  collision  pairs  to  be  considered,  resulting  in  an  over-collision 
ratio  greater  than  one,  as  well  as  a  large  number  of  computationally  expensive  colli¬ 
sions  to  process. 

The  second  option  maintains  a  reasonable  value  of  Np/Nc,  but  results  in  a  very 
small  CEL  number,  causing  an  inefficient  advancement  of  the  solution  and  accumu¬ 
lation  of  statistics.  This  inefficiency  springs  from  the  fact  that  a  small  CEL  number 
implies  that  most  particles  will  require  many  time  steps  to  cross  a  given  cell.  The 
collision  phase  of  each  time  step  will  then  involve  essentially  the  same  group  of  parti¬ 
cles  as  the  previous  time  step  because  very  few  particles  leave  or  enter  the  cell.  This 
situation  is  therefore  equivalent  to  using  a  large  time  step  but  adding  the  compu¬ 
tationally  pointless  exercise  of  moving  the  particles  some  small  distance  at  several 
points  in  the  collision  phase.  It  may  therefore  be  argued  that  it  only  makes  sense  to 
maintain  a  constant  CEL  number,  regardless  of  the  cell  Knudsen  number. 

A  solution  to  this  problem  was  proposed  by  Bartel  et  a/.  [25],  who  recognized  that 
many  of  the  large  number  of  collisions  called  for  by  Equation  4.4  serve  no  purpose 
other  than  to  reinforce  a  Maxwellian  distribution  amongst  the  particles  in  a  given 
cell.  It  should  therefore  be  sufficient  to  restrict  the  number  of  collisions  to  some  small 
value  (comparable  to  the  number  of  particles  in  the  cell)  and  still  take  large  time 
steps  during  the  computation.  In  the  current  terms,  Bartel’s  suggestion  was  to  limit 
the  over-collision  ratio  to  some  (arbitrary)  value  less  than  its  “true”  value  given  by 
Equation  4.4.  This  approach  was  applied  by  Bartel  et  al.  to  a  Couette  flow  and  an 
expanding  nozzle  flow  with  good  results. 


4.2  Numerical  Investigation 

To  explore  these  issues  in  detail,  two  canonical  cases  were  run  with  the  current  code: 
a  Rayleigh  flow  and  a  free  shear  layer.  These  cases  were  chosen  for  their  simplicity, 
the  existence  of  either  analytical  or  published  solutions,  and  their  relatively  small 
CPU  time  requirements.  Ni’s  bump  was  also  run  to  confirm  assertions  made  in  the 
course  of  this  work. 

4.2.1  Rayleigh  Flow 

One-dimensional,  unsteady  Rayleigh  flow  was  chosen  as  the  first  model  problem  for 
this  investigation.  This  flow,  illustrated  in  Figure  4-1,  consists  of  a  stationary  gas 
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subjected  to  the  sudden  acceleration  of  its  lower  boundary  to  a  constant  speed,  Uo- 
It  is  well-suited  to  this  investigation  because  its  length  scales  are  not  imposed  by 
geometry ,  but  rather  by  time  and  viscosity.  In  addition,  it  is  a  one-dimensional 
flow  so  computations  are  compact  and  run  quickly,  allowing  several  test  cases  to  be 
considered. 


Figure  4-1:  Schematic  of  a  One-Dimensional  Rayleigh  Flow 

For  this  series  of  calculations,  the  wall  velocity,  Uoi  was  set  to  0.2  (Mo  =  0.22)  in 
order  to  ensure  the  applicability  of  the  incompressible  Navier-Stokes  solution  (given 
below).  All  cases  were  run  with  Helium  {u  =  0.20)  and  the  CFL  number  was  held  at 
0.3. 

Analytical  Solution 

The  analytical  solution  of  the  Rayleigh  problem  has  two  distinct  regimes  based  on 
the  Knudsen  number.  For  large  Kn,  the  collisionless  Boltzmann  equation  is  applied. 
The  resulting  solution  is  given  by  [9]: 

“  =  (7m-) 

where  R  is  the  gas  constant,  T  is  the  temperature,  and  erfc()  is  the  complimentary 
error  function. 

For  small  Kn,  the  incompressible  Navier-Stokes  equation  is  applied.  For  the 
Rayleigh  problem,  this  reduces  to: 


du  d^u 

dt  dy^ 


(4.7) 


where  u  is  the  kinematic  viscosity  coefficient,  which  is  proportional  to  the  mean  free 
path. 

This  expression  can  be  solved  using  the  slip-flow  boundary  condition  of  Equa- 
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tion  3.1,  resulting  in: 


where  is  a  similarity  variable; 


u  =  Uo 


erfc(^) 
Kn  +  1 


y 

2\/i4' 


(4.8) 

(4.9) 


and  Kn  is  the  Knudsen 
manner  as: 


number,  defined  here  in  a  somewhat 


Kn^ 


X 

s/wnt 


unusual,  but  convenient 


(4.10) 


This  is  a  more  general  version  of  the  classical  Rayleigh  solution. [26]  Note  that, 
in  this  solution.  Kn  is  a  function  of  time  and  becomes  infinite  as  t  — )■  0.  This 
makes  intuitive  sense  because  the  characteristic  scale  of  the  solution  is  the  momentum 
thickness  of  the  viscous  layer,  which  is  initially  zero,  but  grows  with  time.  One 
implication  of  this  time-dependent  Knudsen  number,  however,  is  that  care  must  be 
taken  in  evaluating  the  solution  at  ^  =  0  and  oo,  where  appropriate  limits  of  both  t 
and  y  must  be  computed. 


Well-Resolved  Computation 


This  section  presents  DSMC  results  for  well-resolved  cases  (i.e.  where  Kuc  >  1). 
These  are  intended  to  demonstrate  the  resolution  of  the  current  grid,  which  was  used 
for  all  Rayleigh  flow'  cases,  as  well  as  to  serve  as  points  of  comparison  for  later  sections, 
where  the  geometries  are  not  fully  resolved  and  scaling  issues  are  important. 

The  first  of  these  results,  Figure  4-2,  shows  the  computed  near-wall  velocity  pro¬ 
file,  u{y),  at  t  =  0.004.  Because  this  time  is  significantly  smaller  than  the  mean 
collision  time  (v^/2,  under  the  non-dimensionalizations  of  Section  2.1),  the  collision¬ 
less  Boltzmann  solution  (Equation  4.6)  applies.  The  solid  line  represents  the  analytic 
solution  while  the  symbols  represent  the  DSMC  result.  For  this  computation,  the 
cell  Knudsen  number  was  1.3  x  10"*  and  the  time-step  2.5  x  10“^.  The  agreement  is 
excellent,  providing  evidence  that  the  grid  is  sufficiently  dense  to  accurately  resolve 
the  profile. 

As  the  boundary  layer  grows  with  time,  Kn  increases  through  the  point  where 
collisions  become  important  and  the  fluid  begins  to  behave  as  a  continuum.  When 
Kn  leaves  the  transition  regime  (defined  in  Section  1.2  as  0.1  <  An  <  3),  the  slip- 
flow  Navier-Stokes  solution  of  Equation  4.8  becomes  applicable.  Figure  4-3  contains 
a  comparison  of  DSMC  results  with  this  solution  at  t  =  100.  It  may  be  noted  that 
Kn  =  0.07,  so  there  is  still  a  perceptible  velocity  slip  at  the  wall  (ie.  u(0)  <  Uo). 

If,  as  Figure  4-3  suggests,  DSMC  is  correctly  solving  the  Navier-Stokes  equations, 
the  velocity  distribution  should  be  a  slight  perturbation  from  Maxwellian,  given  by 
the  Chapman-Enskog  distribution  for  a  one-dimensional  isothermal  shear  flow  [1]  as: 


f{u',  v')  =  fo 


uu'v'  du\ 


(4.11) 
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Figure  4-2:  Comparison  of  DSMC  and  collisionless  Boltzmann  solutions  to  the 
Rayleigh  problem  at  t  =  0.004. 


Figure  4-3:  Comparison  of  DSMC  and  slip-flow  Navier-Stokes  solutions  to  the 
Rayleigh  problem  at  t  =  100. 
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where  fo  is  the  Maxwellian  distribution.  This  distribution  was  computed  and  is  shown 
in  figure  4-4,  sampled  at  y  =  2.33,  t  =  40.  Here,  the  mean  velocity,  u,  has  been 
subtracted  and  the  four  quadrants  of  (u',  v')  have  been  collapsed  onto  one  through 
algebraic  operations  which  reinforce  the  perturbation  by  taking  advantage  of  its  anti¬ 
symmetry.  It  should  be  noted  that  the  Chapman-Enskog  distribution  is  formally  only 
valid  for  small  perturbations.  Because  this  is  a  low  Mach  number  computation,  the 
thermal  fluctuations  are  significant  and  therefore  only  qualitative  comparisons  are 
appropriate.  Nevertheless,  the  deviation  from  Maxwellian  computed  in  the  DSMC 
simulation  is  in  good  agreement  with  its  predicted  structure,  indicating,  as  expected, 
that  the  gas  is  weakly  perturbed  from  equilibrium  by  the  shear. 


Figure  4-4:  Distribution  of  velocity  perturbations  for  a  well-resolved  computation. 

In  the  following  section,  where  this  case  is  recomputed  with  various  values  of  Kric, 
a  quantitative  comparison  of  the  shear  layers  will  be  needed.  One  such  a  measure 
is  obtained  by  performing  a  non-linear  least-squares  fit  of  the  streamwise  velocity 
data  to  the  profile  of  Equation  4.8,  with  ly  as  the  free  parameter.  To  demonstrate  the 
dependability  of  this  measure.  Figure  4-5  shows  the  time-dependence  of  the  computed 
viscosity  for  the  well-resolved  computation  as  it  evolves  from  t  =  0  to  200. 

With  the  non-dimensionalizations  used  in  the  code,  the  normalized  kinematic 
viscosity  coefficient  of  Helium  is  0.64.  From  Figure  4-5,  it  may  be  seen  that  the  DSMC 
solution  under-predicts  this  viscosity  for  small  times,  but  appears  to  asymptote  to 
nearly  the  correct  value  as  Kn  — >  0.  The  under-prediction  at  early  times  may  be 
explained  by  noting  that  Kn  is  not  yet  in  the  region  of  validity  for  Equation  4.8.  The 
slight  over-prediction  at  later  times  is  most  likely  due  to  a  variety  of  factors  including 
statistical  scatter  due  to  the  low-Mach  number  and,  a  weak  influence  of  the  far  wall 
(located  200  A  from  the  moving  surface). 
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Figure  4-5:  DSMC-computed  viscosity  and  Knudsen  number  as  a  function  of  time  for 
a  well-resolved  computation. 

Under- Resolved  Computations 

The  previous  section  demonstrated  that  the  DSMC  code  can  accurately  reproduce  the 
analytical  solution  of  the  Rayleigh  problem.  This  is  not  surprising,  however,  because 
the  method’s  ability  to  model  the  true  gas  physics  was  demonstrated  in  the  previous 
chapter  and  by  many  other  researchers.  These  results  were  presented  mainly  as  a 
means  of  verifying  the  proper  matching  of  the  code  and  grid  to  the  problem  and  for 
calibrating  the  measures  used  in  the  scaling  investigation. 

This  investigation  is  begun  in  the  current  section  by  manipulating  the  previous 
Rayleigh  calculation.  In  all  cases,  the  grid  geometry  was  maintained,  but  it’s  lin¬ 
ear  dimensions  were  increased  by  a  multiplicative  factor.  To  hold  the  CFL  number 
constant^,  this  factor  was  also  applied  to  the  time  step.  The  number  of  particles  was 
not  changed. 

In  a  scaled  calculation  with  a  high  over-collision  ratio,  it  is  expected,  as  argued 
above,  that  the  particles  in  each  cell  will  approach  a  Maxwellian  distribution,  rather 
than  the  proper,  perturbed  Maxwellian  shown  in  Figure  4-4.  This  hypothesis  was 
tested  by  running  a  case  with  an  over-collision  ratio  of  approximately  5  {Kric  =  0.05). 
The  distribution  was  then  sampled  from  the  same  grid  location  as  the  previous  case 
at  a  time  selected  to  obtain  approximately  the  same  point  in  the  profile  (ie.  the 
same  value  of  u{y)/Uo).  With  the  current  scaling,  this  corresponds  to  y  =  46.67 
and  t  =  2700.  The  resulting  distribution  is  shown  in  Figure  4-6,  collapsed  onto  one 


^Selected  computations  performed  with  lower  CFL  numbers  yielded  almost  identical  results, 
supporting  the  assertions  made  in  Section  4.1. 
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quadrant  in  the  same  manner  used  in  Figure  4-4. 


Figure  4-6:  Distribution  of  velocity  perturbations  for  an  over-collided  (under-resolved) 
computation 

In  contrast  to  the  previous  case,  there  is  no  well-defined  structure  to  the  pertur¬ 
bations  in  Figure  4-6.  This  lack  of  structure  may  be  quantified  by  noting  that  the 
integral  of  the  collapsed  distribution  shown  here  is  an  order  of  magnitude  smaller 
than  that  for  the  well-resolved  case.  The  implications  of  this  result  are  somewhat 
subtle:  if  each  cell  in  the  over-collided  computation  contains  an  equilibrium  distri¬ 
bution,  then  the  Navier-Stokes  equations,  which  correspond  to  a  weak  perturbation 
from  the  Maxwellian  state,  are  not  being  solved,  but  rather  the  Euler  equations,  which 
represent  an  ideal  gas  in  a  perpetual  state  of  equilibrium.  Indeed,  the  closer  each  cell 
is  driven  to  equilibrium,  the  more  “ideal”  the  fluid  should  become.  If  this  is  in  fact 
true,  one  could  argue,  then  the  computed  viscosity  of  the  Rayleigh  layer  should  go  to 
zero  as  the  cell  Knudsen  number  decreases,  confining  the  effect  of  the  wall  motion  to 
an  increasingly  narrow  layer  which  asymptotes  to  a  vortex-sheet  singularity  at  ?/  =  0. 

When  this  assertion  is  tested  through  a  series  of  runs  with  different  grid  scalings, 
however,  a  surprising  result  is  found:  instead  of  a  decrease  in  viscosity  as  Kric  ->  0, 
the  opposite  occurs.  As  Figure  4-7  indicates,  the  gas  viscosity  calculated  from  the 
DSMC  computations  increases  as  a  function  of  the  cell  size  (Aa;  =  IjKn^. 

Three  sets  of  data  are  reported  in  Figure  4-7:  in  the  first  set,  denoted  by  circles, 
the  number  of  collisions  computed  during  each  time  step  was  not  limited;  in  the  second 
set,  denoted  by  stars,  a  collision  limiter  of  5  was  enforced;  in  the  last  set,  denoted  by 
crosses,  a  collision  limiter  of  1  was  employed.  (A  collision  limiter  of  1  implies  that 
A/p  =  Ac  regardless  of  the  value  of  called  for  by  the  NTC  equation  (Equation  4.4).) 
At  low  values  of  1/Ahc,  the  correct  v  is  produced.  However,  as  the  cell  size  rises 
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Figure  4-7:  DSMC-computed  viscosity  versus  cell  size. 

above  about  10,  the  eflfective  viscosity  begins  to  increase.  This  increase  appears  to 
be  linear  over  a  very  wide  range  of  Kric- 

This  behavior  stems  from  the  particulate  nature  of  the  simulation.  At  low  Mach 
numbers,  the  random  thermal  velocity  of  the  gas  will  always  result  in  particles  mov¬ 
ing  in  the  y-direction  even  though  there  is  no  net  motion  normal  to  the  surface.  In 
addition,  because  DSMC  chooses  collision  partners  in  a  given  cell  without  regard  to 
their  location,  momentum  is  uniformly  diffused  through  the  cell  within  a  few  time- 
steps.  In  a  properly-resolved  computation,  where  the  cell  dimension  is  approximately 
one  mean-free-path,  this  behavior  is  molecularly  “correct”  and  results  in  the  proper 
physical  viscosity  of  the  fluid  (helped,  of  course,  by  an  appropriate  value  of  the  VHS 
exponent,  u).  However,  in  an  under-resolved  computation,  the  disregard  for  a  par¬ 
ticle’s  location  in  the  cell  results  in  a  diffusion  of  momentum  due  to  thermal  motion 
that  is  far  in  excess  of  the  physical  viscosity  and  it  is  this  artificial  viscosity  that  is 
exhibited  in  these  results.  In  many  cases,  this  thermal  motion  would  not  present  a 
problem.  However  in  the  presence  of  a  strong  mean  velocity  gradient  (as  in  the  case 
of  the  wall-driven  Rayleigh  problem),  the  thermal  motion,  coupled  with  the  mean 
shear,  results  in  an  artificial  Newtonian  viscosity,  inversely  proportional  to  the  cell 
Knudsen  number,  Knc,  and  greater  than  the  physical  viscosity. 

This  argument  suggests  that  the  artificial  viscosity  should  be  most  apparent  at 
low  Mach  number  and  will  abate  as  Ma  increases  because  the  relative  importance  of 
thermal  fluctuations  is  reduced.  In  other  words,  at  high  speeds,  particles  will  be 
swept  from  the  domain  by  the  mean  flow  before  they  have  the  opportunity  to  diffuse 
appreciably  in  the  transverse  direction.  This  was  not  found  to  be  the  case  for  the 


50 


Rayleigh  problem,  however.  The  apparent  discrepancy  is  explained  by  noting  the 
flow’s  one-dimensional  nature,  which  was  exploited  in  this  calculation  by  employing 
a  Svrap-around’  boundary  condition.  Here,  exiting  particles  are  reintroduced  to  the 
domain  on  the  opposite  edge,  where  they  can  continue  to  diffuse  in  the  y-direction. 
A  2-D,  spatially-developing  flow  is  therefore  needed  to  verify  the  above  argument. 


4.2.2  Boundary- Free  Shear  Flows 

A  boundary-free  shear  flow  was  selected  to  test  the  above  assertion.  Here,  two  streams 
at  different  velocities  are  introduced  at  x  =  0;  the  upper  stream  moving  at  and 
the  lower  at  U2-  This  arrangement  is  illustrated  in  Figure  4-8. 


Figure  4-9  shows  contours  of  velocity  in  the  x  —  y  plane  for  a  well-resolved  calcu¬ 
lation  where  U\  —  1.5  (Ma  =  1.6),  U2  =  2.0  (Ma  =  2.2),  and  the  working  fluid  is 
Argon  (tu  =  0.31).  Note  that  the  singularity  at  x  =  0  quickly  diffuses  and  a  slowly 
thickening  shear  layer  develops. 

By  integrating  the  streamwise  velocity  across  the  layer,  a  flow  can  be  characterized 
by  its  momentum  thickness,  d\ 

9  =  -n)(u-U2)dy  (4.12) 

where  9  has  been  normalized  by  the  upper  stream  velocity,  Ui  and  the  computational 
domain  height,  H.  This  momentum  thickness  is  shown  in  Figure  4-10  as  a  function 
of  X  for  several  computations  at  different  Kric- 

Theoretically,  the  spatial  evolution  of  a  laminar  shear  layer  shows  a  square-root 
growth  in  momentum  thickness  [26]: 


9  oc  i/x 


(4.13) 
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Figure  4-9:  Contours  of  streamwise  velocity  in  a  well-resolved  DSMC  computation  of 
a  supersonic  shear  flow. 

which  is  well  captured  by  all  of  the  test  cases.  Anamolies  are  observed,  however, 
when  comparing  these  cases  to  one  another. 

First,  comparing  the  well-resolved  case  (120x60)  to  the  first  scaled  case  (2400x 
1200),  it  is  expected  that  the  twenty- fold  increase  in  linear  dimensions  will  result 
in  a  \/^  =  4.47  decrease  in  the  shear  layer  thickness  at  the  same  value  of  x/H. 
In  Figure  4-10,  however,  the  shear  layer  only  decreases  by  a  factor  of  approximately 
1.3.  This  is  again  explained  by  the  artificial  viscosity  introduced  by  the  under-resolved 
grid:  as  the  grid  is  scaled  from  the  well-resolved  case,  this  artificial  viscosity  surpasses 
the  physical  viscosity,  altering  the  proper  relationship  between  the  curves  because  the 
cases  appear  to  be  run  with  different  fluids. 

Second,  a  further  doubling  from  a  scale  factor  of  20  to  40  (relative  to  the  base¬ 
line  case)  is  expected  to  further  reduce  the  shear  layer  thickness,  but  was  found 
to  have  a  negligible  influence.  This  “saturation”  effect  is  related  to  the  particle- 
transport/  collisional-smearing  mechanism  causing  the  artificial  viscosity.  Because  the 
grid  dimensions  and  time  step  were  scaled  by  the  same  factor  (ie.  the  CFL  number 
was  held  constant)  identical  particles  moving  through  the  two  cases  will  pass  through 
the  same  locations  on  the  grid  and  undergo  the  same  number  of  collision  phases.  In 
other  words,  for  this  Mach  number,  a  characteristic  level  of  artificial  viscosity  has 
been  reached  which  is  dependent  on  the  grid  itself.  For  these  scaling  factors,  this 
grid-based  viscosity  is  so  much  greater  than  the  physical  viscosity  that  the  flow  has 
become  insensitive  to  its  physical  dimensions. 

Finally,  when  the  Mach  number  of  the  incoming  flow  is  increased,  keeping  the 
ratio  U1/U2  constant,  the  theory  predicts  no  change  in  the  shear  layer  thickness.  The 
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Figure  4-10:  Development  of  shear  layer  momentum  thickness  as  a  function  of  down¬ 
stream  distance. 

trace  in  Figure  4-10  for  this  case  (labelled  ‘4800x2400  High  Speed)  shows  a  thinner 
shear  layer,  however  (still  not  reaching  its  correct  value,  but  closer  nonetheless).  This 
improvement  is  due  to  the  lower  importance  of  thermal  velocities  with  respect  to  the 
stream  velocity,  which  results  in  a  lower  artificial  viscosity.  It  should  be  noted  that, 
as  in  the  Rayleigh  problem  presented  in  the  previous  section,  the  artificial  viscosity 
was  observed  to  be  relatively  insensitive  to  the  use  of  a  collision  limiter  and  to  the 
choice  of  the  CFL  number.  In  this  geometry,  however,  it  was  observed  that  a  collision 
limiter  did  thin  the  layer  very  slightly.  This  is  perhaps  due  to  the  fact  that,  with  a 
collision  limiter,  there  is  slightly  less  “uniformity”  within  each  cell  because  there  are 
fewer  collisions,  so  momentum  is  not  smeared  across  the  shear  layer  as  rapidly. 

4.2.3  Euler  Flow 

The  results  of  the  preceding  section  illustrate  two  competing  trends  in  DSMC  com¬ 
putations  as  the  cell  Knudsen  number  decreases.  On  one  hand,  a  high  value  of  the 
over-collision  ratio,  Np/Nc,  drives  the  particles  in  each  cell  toward  a  Maxwellian,  or 
equilibrium,  distribution.  If  this  were  the  only  effect  present,  DSMC  would  therefore, 
in  these  cases,  formally  solve  the  Euler  equations,  which  are  the  governing  dynamic 
equations  for  an  equilibrium  gas.  A  competing  effect  exists,  however,  because  the 
essential  nature  of  DSMC  -  the  random  motion  of  particles  coupled  with  a  disregard 
for  the  location  of  collision  partners  in  a  cell  -  creates  an  artificial  viscosity  from  the 
diffusion  of  momentum  through  the  domain. 

From  this,  one  might  expect  that  a  scaled  DSMC  code  should  fail  most  dramati- 
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cally  when  attempting  a  low  Mach  number  viscous  flow  and  be  optimal  for  modeling 
a  high  Mach  number,  inviscid  flow.  The  former  assertion  was  demonstrated  in  the 
previous  .sections;  this  section  presents  results  to  confirm  the  latter. 

The  test  case  selected  for  this  task,  known  as  “Ni’s  Bump”,  is  often  run  to  demon¬ 
strate  inviscid  CFD  methods.  [27]  The  geometry  consists  of  a  two-dimensional  duct 
with  a  circular-arc  excrescence  placed  on  the  floor.  A  bump  with  20%  of  the  duct 
height  was  chosen  for  this  work.  This  geometry  and  a  traditional  Euler  code  solution 
are  shown  in  Figure  4-11.  The  inlet  Mach  number  is  3.28  and  the  working  fluid  is 
Argon. 


Results 

Figure  4-12  shows  Mach  Number  contours  for  a  well-resolved  DSMC  computation  of 
Ni’s  bump  with  duct  dimensions  of  125x30  A.  This  run  utilized  about  42,000  particles 
on  a  grid  of  1,860  cells  with  u)  =  0.31.  The  flow  reached  statistical  equilibrium 
in  approximately  10  minutes,  but  an  additional  two  hours  were  required  to  gather 
adequate  statistics  for  data  analysis  (on  an  SGI  Indigo).  Although  specular  walls 
were  employed,  and  thus  a  viscous  solution  is  not  being  attempted,  it  is  instructive 
to  note  that  the  effective  Reynolds  number  of  the  DSMC  computation,  based  on  the 
inlet  Mach  number  and  the  bump  length  (which  is  dimensional  due  to  the  mean  free 
path  grid  scaling),  is  243. 


Bump. 

By  comparing  this  solution  to  that  of  the  Euler  code,  it  may  be  concluded  that 
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Figure  4-13:  Mach  number  contours  of  a  scaled  DSMC  computation  of  Ni’s  bump. 


the  flow  is  faithfully  reproduced  except  for  an  apparently  large  shock  thickness  in 
the  DSMC  result.  This  is  to  be  expected  since,  although  there  are  no  viscous  effects 
due  to  the  boundaries,  the  DSMC  method  is  accurately  solving  the  shock  structure, 
where  viscous  effects  are  strong.  Because  the  entire  grid  is  only  125  mean  free  paths  in 
extent,  however,  this  shock  appears  excessively  thick,  though  its  thickness  is  actually 
quite  comparable  to  the  physical  case  («  10  A).  In  addition,  other  features  of  the  flow, 
including  the  shock  angles,  the  leading  and  trailing  edge  shocks,  the  shock  reflections 
and  the  shock-shock  interaction  toward  the  rear  of  the  duct,  are  properly  reproduced. 

Figure  4-13  shows  the  same  computation,  scaled  by  a  factor  of  1,000  so  the  di¬ 
mensions  of  the  duct  are  now  125,  000  x  30,  000  A  (approximately  9x2  millimeters  at 
atmospheric  conditions).  For  this  case,  a  collision  limiter  of  one  was  used.  The  grid 
resolution  was  also  doubled  in  each  direction  (compared  to  the  last  case)  to  7,750 
cells.  This  was  necessary  because  the  shock  thickness  was  observed  to  asymptote  to  a 
larger-than-expected  minimum  value  as  the  cell  Knudsen  number  decreased.  This  was 
identified  as  a  grid  resolution  issue  and  it  was  found  that,  in  common  with  standard 
CFD  practice,  the  computed  flow  converged  to  a  grid-independent  solution  through 
successive  refinements  in  the  grid  (the  average  number  of  particles  per  cell  was  held 
constant). 

When  compared  with  the  previous  case,  the  shock  appears  to  be  much  thinner, 
although  it  is  now  actually  thicker  than  a  physical  shock  due  to  the  much  larger  grid 
dimensions.  This  large  shock  thickness  is  due  to  the  artificial  viscosity  inherent  in  a 
continuum  application  of  DSMC.  This  does  not  render  the  solution  invalid,  however;  it 
is,  in  fact,  common  practice  in  inviscid  continuum  CFD  codes  to  explicitly  introduce 
an  artificial  viscosity  (usually  second-order)  in  order  to  stabilize  the  solution  near 
steep  gradients  such  as  shocks.  The  effect  of  this  artificial  viscosity  is,  as  in  the 
DSMC  solution,  a  thickening  of  the  shocks.  DSMC  does  not  require  any  explicit 
introduction  of  artificial  viscosity,  however;  it  provides  its  own,  as  demonstrated  by 
these  results.  Additionally,  no  oscillations  in  the  solution  are  observed  in  the  regions 
surrounding  the  shock. 
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Chapter  5 
Conclusions 


This  report  has  presented  many  aspects  of  DSMC’s  application  to  micromechanical 
devices.  In  the  preceding  chapters,  the  motivation  behind  this  application,  the  al¬ 
gorithm  developed  for  it,  some  of  its  results,  and  the  issues  involved  with  scaling  it 
beyond  its  traditional  range  were  discussed.  This  chapter  draws  conclusions  from 
these  efforts  to  evaluate  DSMC’s  potential  as  a  tool  for  MEMS  development. 

As  discussed  in  Chapter  1,  the  high  Knudsen  numbers  of  many  flows  of  interest 
to  MEMS  designers  make  them  inaccessible  to  continuum-based  numerical  methods. 
While  analytical  models  are  being  developed  for  moderately-rarefied  flows  in  simple 
geometries  [23],  the  complexity  and  mixed-regime  nature  of  many  MEMS  limits  the 
applicability  of  these  methods.  Experimental  study  of  these  flows  has  also  proven  to 
be  difficult  due  to  the  miniscule  quantities  which  must  be  measured.  In  the  course 
of  their  micro-channel  investigation,  for  example,  Arkilic  et  al.  [23]  were  required  to 
develop  a  system  capable  of  measuring  mass  flows  on  the  order  of  nanograms  per 
second.  In  response  to  this  situation.  Chapter  3  presented  a  small  subset  of  the 
micro-flows  of  interest,  to  both  the  fluid-dynamicist  and  the  MEMS  designer,  which 
are  easily  treated  with  DSMC. 

In  Section  3.1,  DSMC  results  for  a  slip-flow  regime  micro-channel  were  compared 
to  an  analytical  solution  presented  by  Arkilic  et  a/.  [23]  Several  aspects  of  the  computed 
and  theoretical  streamwise  velocity  and  pressure  distributions  were  considered.  In 
all  cases,  the  numerical  results  were  found  to  compare  well  with  their  analytical 
predictions.  In  addition,  by  manipulating  the  analytical  relations,  an  expression  for  a 
streamwise  velocity  ‘similarity’  profile  was  developed.  It  was  subsequently  found  that 
the  computed  velocity  distribution,  which  is  a  function  of  x,  collapsed  quite  well  to 
the  predicted  x-independent  profile  when  the  flow  variables  were  processed  according 
to  the  theoretical  expression. 

These  findings  support  the  accuracy  of  both  the  code  and  the  analytical  work 
and  demonstrate  another  valuable  function  of  DSMC  computations:  the  evaluation 
of  analytical  models  for  rarefied  flow  behavior.  Developing  these  models  is  an  impor¬ 
tant  task  because  a  great  number  of  MEMS  geometries  fall  into  the  slip-flow  regime. 
Analytical  tools  have  great  value  to  the  designer  because  they  are  much  less  expensive 
than  a  full  DSMC  calculation  but  are  able,  for  certain  geometries  and  flow  regimes, 
to  quantify  behavior  invisible  to  continuum  techniques.  They  are  difficult  to  validate, 
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however,  because,  as  mentioned  above,  the  effects  they  describe  are  often  too  small 
for  effective  experimental  investigation. 

In  Section  3.2,  a  DSMC  solution  for  a  transition  regime  micro-channel  was  com¬ 
pared  to  predictions  from  the  slip-flow  expressions  of  Arkilic  et  al.  By  considering 
the  form  of  disagreement  between  these  results,  several  aspects  of  flow  behavior  in 
this  regime  were  illuminated.  This  is  an  especially  important  application  of  DSMC 
because  very  few  techniques,  either  analytical  or  numerical,  are  available  for  analyz¬ 
ing  these  flows.  This  is,  consequently,  the  least  understood  of  the  four  /sfn-regimes, 
posing  a  significant  problem  for  the  MEMS  community  because  the  geometries  and 
working  conditions  for  many  devices  place  them  here.  Additionally,  through  the  sim¬ 
ilarity  analysis  developed  in  the  previous  section,  the  onset  of  transition  behavior  was 
observed  as  the  Knudsen  number  increased  down  the  channel.  This  type  of  work  is 
intended  to  determine  the  applicable  range  of  various  models  and  the  mode  of  their 
failure.  This  is  valuable  because  applying  convenient  models  to  the  largest  possible 
range  of  problems  is  preferable  to  performing  large  numerical  simulations,  especially 
early  in  the  design  process. 

In  the  final  section  of  Chapter  3,  parabolic  micro-nozzles  were  examined.  These 
cases  are  intended  to  represent  devices  which  contain  complexities  not  amenable  to 
other  forms  of  solution.  Their  more  complicated  geometry,  lack  of  isothermal  flow, 
and  sharp  gradients  pose  serious  problems  for  many  types  of  analysis.  No  special 
considerations  for  these  features  were  required  to  perform  the  DSMC  calculation, 
however.  This  demonstrates  the  versatility  of  the  method  and  its  value  for  complex 
flows,  which  are  common  in  MEMS  devices  and  will  become  more  so  as  the  field 
matures. 

Chapter  4  extended  this  work  by  investigating  the  issues  inherent  in  scaling  a 
DSMC  code  to  treat  continuum  flows.  It  is  important  to  understand  these  issues 
because  many  MEMS  devices  contain  regions  in  this  regime.  Due  to  its  molecular- 
level  scaling,  however,  treating  them  with  a  fully-resolved  DSMC  calculation  is  very 
computationally  expensive.  A  greater  understanding  of  DSMC’s  behavior  when  its 
scaling  constraints  are  relaxed  is  therefore  valuable  because  it  may  enable  the  simpli¬ 
fied  treatment  of  certain  geometries. 

In  this  chapter,  two  primary  effects  were  observed  when  the  grid  was  scaled  beyond 
the  traditional  constraints  of  DSMC  calculations.  The  first  of  these  was  previously 
identified  by  Bartel  et  al.  [25]:  when  the  scaling  reaches  a  point  where  particles 
undergo  multiple  collisions  every  time  step,  the  velocity  distribution  in  each  cell 
approaches  an  equilibrium  (Maxwellian)  form.  The  implication  of  this,  which  was 
not  previously  realized,  is  that  the  DSMC  technique  is  formally  solving  the  Euler 
equations,  not  the  Navier-Stokes  equations,  in  these  regions.  In  cases  where  this  is 
appropriate,  Bartel  et  al.  note  that  the  computation  may  be  greatly  accelerated  by 
employing  a  ‘collision  limiter’  large  enough  to  enable  each  cell  to  reach  equilibrium, 
but  not  so  large  that  this  state  is  needlessly  reinforced  through  continued  interaction. 

A  second  implication  of  DSMC’s  application  to  continuum  flows  is  that  the  thermal 
velocity  of  the  computational  particles,  coupled  with  the  physically  large  cell  size, 
results  in  an  artificial  viscosity  which  is  proportional  to  the  cell  Knudsen  number 
and  decreases  with  flow  Mach  number  (because  the  thermal  fluctuations  become  less 
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important).  This  is  a  consequence  of  the  DSMC  collision  algorithm,  which  tends  to 
distribute  momentum  and  energy  uniformly  through  the  entire  cell  due  to  its  disregard 
for  the  position  of  colliding  partners  within  that  cell.  This  tendency  is  particularly 
strong  when  a  large  number  of  collisions  are  performed  during  each  time  step,  as  is 
the  case  when  Kric  is  small.  This  artificial  viscosity  can  be  a  nuisance  when  solving 
the  Navier-Stokes  equations  at  low  Reynolds  numbers,  where  the  physical  viscosity 
is  important  to  the  overall  flow  behavior,  because  it  is  implicit  to  the  scheme  and 
represents  a  first-order  error  term.  Large  reductions  must  therefore  be  made  in  the 
cell  size  to  bring  the  computed  viscosity  acceptably  close  to  the  true  value.  This 
artificial  viscosity  may  be  an  asset  for  solving  Euler  flows,  however,  because  it  acts 
as  a  built-in,  self-adjusting  means  of  broadening  poorly-resolved  flow  features  such 
as  shock  structures,  serving  a  similar  function  as  the  artificial  viscosity  explicitly 
included  in  many  traditional  CFD  methods. 

In  summary,  the  unique  features  of  rarefied  flows  have  many  implications  for 
MEMS  designers.  The  larger  mass  flow  rate  for  a  given  geometry  and  inlet  conditions, 
found  in  Sections  3.1  and  3.2,  must  be  considered  when  designing  control  systems  for 
micro-chemical  reactors  and  the  decreased  thermal  communication  between  the  fluid 
and  its  boundaries,  demonstrated  in  Section  3.3,  must  be  considered  in  applications 
involving  temperature  measuring  devices  or  heaters,  for  example.  These  flows  are 
not  well-studied,  however,  due  to  the  breakdown  of  the  transport  assumptions  which 
characterize  the  Navier-Stokes  equations,  the  basis  of  most  common  fluid-dynamic 
analysis  tools. 

DSMC’s  ability  to  calculate  in  any  of  the  four  iifn-regimes  without  modification 
makes  it  ideal  for  investigating  flows  related  to  MEMS  devices,  particularly  those  with 
regions  in  different  regimes.  It  is  therefore  useful  to  MEMS  researchers  for  validating 
analytical  models,  investigating  flow  behavior,  and  modeling  potential  devices.  In 
addition,  it  is  quite  straightforward  to  include  flow  features  such  as  chemical  reactions, 
multi-species  mixing  and  particle  transport  in  the  code  due  to  its  particulate  nature 
and  modular  construction.  In  addition,  inclusion  of  unstructured  grid  capability 
and  the  trajectory-tracing  movement  scheme  enable  the  code  to  handle  arbitrary 
geometries.  This  is  a  valuable  asset  when  analyzing  the  complex  structures  included 
in  many  of  these  devices.  Finally,  its  cell-based  structure  makes  it  well-suited  for 
parallel  computation,  which  is  an  increasingly  important  attribute  for  a  large-scale 
numerical  scheme  on  modern  computers. 

These  attributes  are  especially  important  to  consider  in  light  of  the  scaling  inves¬ 
tigation  performed  in  Chapter  4.  Given  the  speed  and  maturity  of  continuum  CFD 
methods,  it  appears,  at  first  glance,  that  the  scaling  attempted  in  this  chapter  is 
irrelevant  because,  given  the  choice,  the  continuum  method  would  always  be  selected. 
This  is  not  necessarily  true,  however,  because  continuum  methods  are  often  more  dif¬ 
ficult  to  parallelize  and  complexities,  such  as  those  mentioned  above,  each  represent 
large  increases  in  computational  work.  It  is  therefore  possible  that  the  speed  of  a 
DSMC  calculation  in  certain  problems  and  on  certain  machines,  may  be  quite  com¬ 
parable  to  that  of  a  continuum  technique.  In  addition,  mixed  regime  flows  are  not 
generally  treatable  with  continuum  methods  because  they  cannot  effectively  model 
the  rarefied  areas.  While  hybrid  continuum/particulate  codes  are  under  development 
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[28]  for  these  flows,  it  would  be  desirable  to  treat  the  entire  domain  uniformly  to  max¬ 
imize  code  flexibility  and  avoid,  if  possible,  the  problems  associated  with  interfacing 
two  techniques  of  entirely  different  natures. 

Many  aspects  of  DSMC’s  application  to  MEMS  remain  to  be  explored  and  devel¬ 
oped.  To  model  actual  devices  in  their  entirety,  for  example,  the  current  code  will 
require  an  upgrade  to  three-dimensional  calculations.  The  corresponding  increase  in 
computational  requirements  will  then  make  it  necessary  to  move  from  the  current, 
workstation-class,  computers  to  large-scale  supercomputers.  As  mentioned  previously, 
DSMC’s  structure  is  very  amenable  to  efficient  parallelization,  so  this  move  will  likely 
be  to  a  parallel  machine.  On  these  architectures,  many  computational  issues  remain 
to  be  explored  which  promise  to  increase  the  speed  and  efficiency  of  the  technique. 
In  complex  cases  for  which  a  particle  description  is  an  asset,  such  as  multispecies  and 
chemically  reacting  flows,  this  method  may  then  be  advantageous  even  in  regimes  to 
which  continuum  techniques  are  currently  applied.  For  mixed-regime  flows,  hybrid 
DSMC-continuum  techniques  represent  a  promising  area  which  still  requires  signif¬ 
icant  investigation.  When  mature,  these  techniques  will  allow  both  methods  to  be 
applied  where  they  are  most  desirable  and  ‘swapped  out’  where  they  become  ineffi¬ 
cient  or  inaccurate.  A  numerical  solution  will  then  involve  applying  the  best  of  each 
discipline  to  a  problem,  rather  than  choosing  one  or  the  other  when  setting  up  the 
calculation  and  patching  together  corrections  and  modifications  for  certain  portions 
of  the  flow. 
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Appendix  A 

Particle  Movement  Function 


This  appendix  contains  the  function  which  performs  particle  movement.  It  is  intended 
to  serve  only  as  an  annotated  listing;  the  full  description  of  the  movement  operation 
is  given  in  Section  2.4. 

Arguments 

The  movement  function  takes  four  arguments.  The  first  argument,  ‘cel’  is  a  pointer 
to  the  current  cell  of  the  particle(s)  to  be  moved.  The  second  argument,  ‘particles’ 
is  a  pointer  to  the  first  particle  in  a  list  of  those  to  be  moved.  In  the  case  of  a 
routine  global  move,  such  as  at  the  beginning  of  a  time  step,  this  is  simply  the 
cell’s  particle  list.  For  single  particle  moves,  such  as  when  one  is  introduced  at  an 
inflow/outflow  face  or  emitted  from  a  diffusely- reflecting  surface,  this  is  simply  the 
address  of  the  particle  to  be  displaced.  The  third  argument,  ‘nparts’  is  simply  the 
number  of  particles  to  be  moved.  The  final  argument  ‘sel_move_flag’,  is  an  integer 
constant  that  tells  the  function  whether  or  not  to  read  the  movement  times  remaining 
for  individual  particles.  When  the  flag  is  set  to  0,  all  particles  are  considered  to  be 
moving  an  entire  time  step,  as  is  the  case  in  a  global  move. 

Cell  Face/Node  Numbering  Scheme 

When  calculating  intersection  times,  a  simple  relation  between  the  face  and  node 
indices  is  required.  In  order  to  preserve  the  trajectory-tracing  method’s  ability  to 
function  in  cells  with  an  arbitrary  number  of  faces,  the  following  scheme  was  adopted: 
face  i  connects  nodes  i  and  {i  +  l)%ns  where  i  is  an  index  which  starts  at  zero,  n*  is 
the  total  number  of  sides  for  the  cell,  and  %  is  the  modulus  operator,  which  returns 
the  remainder  resulting  from  the  division  of  its  operands.  The  macro,  ‘IRING’,  that 
appears  in  this  function  performs  this  modulo-division  operation  to  cycle  through  cell 
nodes.  For  the  current,  triangulated  mesh,  it  is  given  by: 

#define  IRING  (A)  =  ((A)  %  3)) 

This  numbering  scheme  is  shown  in  Figure  A-1. 
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Code  Listing 

void  Move_Particle(struct  cell_unit  *cel.  struct  atom  *particles, 
unsigned  short  n_parts,  short  sel_move_flag) 

short  fc,  face,  entry_face,  yjntersect; 
float  timejnterval,  tjnt,  t_min,  traj_slope; 
unsigned  curr_cell; 
struct  atom  *part; 

if(!sel_move_flag)  { 

j****  Set  particle  movement  time  interval  to  be  the  time  step  ^^**1 
timejnterval  =  T_STEP; 

!****  Set  the  entry  face  to  be  a  nonexistant  side  ****/  lo 

entry_face  =  3; 

} 

for  (part  =  particles;  part  <  particles  +  n_parts;  part+-{-){ 

j  Skip  if  this  is  an  empty  particle  position  ****! 
if  (part“>dest_celi  — =  0)  continue; 

j  ****  If  particle  has  somehow  left  grid,  bring  it  back  ****! 
if  (part->x  >  x_max+0.1  ||  part— >y  >  y_max+0.1  ||  part->x  <-0.1 
II  part->y<-0.1){ 

printf('’Relocating  Particle  at  t  =  '/.f :  \n",t^lobal); 

printf(''  x:  7,f  y:  7,f  u:  7,  .24f  v:  7# .  24f  cell:  7.d  particle:  7.<i\n",  20 
part->x,  part->y,  part->u,  part->v,  cel  -  cell,  part  -  particles); 
Random_Plcmt2(cel,  part); 

} 

if(sel_move_flag){ 

j****  Determine  time  interval  over  which  to  attempt  to  move  particle 
timejnterval  =  part— >mvmt_time; 

I  Skip  if  this  particle  is  finished  moving  ****/ 
if  (timejnterval  ==  0.)  continue; 

} 

/  ****  Determine  which  face  is  ineligible  for  intersection  because  it  was  30 

just  crossed  (selective  moves  only)  ****/ 
if(sel_move_flag) 


62 


entryjace  —  (short )  part— >crosd_face; 

/  Calculate  the  slope  of  the  particle  trajectory 
if(part->u  =—  0.) 

traj.slope  =  INFINITY; 

else 

traj_slope  =  part->v  /  part->u; 
if  (traj_slope  *  traj_slope  >  1.) 

yjntersect  =  1;  40 

else 

yjntersect  =  0; 

f****  Initialize  tjnin  with  a  large  value  ****/ 
t_min  =  1.05  *  T.STEP; 

for  (fc  =  0;  fc  <  3;  fc++){ 

if  (fc  ==  entry_face)  continue; 

If  particle  path  is  parallel  to  face,  they  will  never  intersect  ****! 
if(cel->slope[fc]  ==  traj_slope)  continue; 

!  ****  If  particle  ^s  vertical  component  is  greater  than  its  horizontal  (in 

absolute  value),  use  y~intersects  for  time  determination  ****/  so 

if(yjntersect) 

tjnt  =  (((traj_slope  *cel->slope[fc]  *(cel->x[fc]  -  part->x) 

+  cei->slope[fc]  *part->y  ~traj_slope  *cel->y[fc]) 

/  (cel->slope[fc]  -traj_slope))  -  part->y)  /part->v; 

else 

tJnt  =  (((cel->y[fc]  -  part->y  +  traj_slope*part->x 

-  cel->slope[fc]*cel— >x[fc])  /  (traj_slope  -  cel->slope[fc])) 

—  part— >x)  /  part— >u; 

I  ****  If  t_int<0  then  particle  is  heading  the  wrong  way  to  hit  this  face  ****/ 

if  (tjnt  <  —  I_TEST)  continue;  eo 

I  ****  If  the  intersect  time  is  zero,  make  sure  this  particle  is  heading 
in  a  direction  which  makes  striking  this  face  possible  before  accepting 
it  to  avoid  an  infinite  loop  of  zero  time  steps 
if  (tjnt  <=  I_TEST  &&  ((part->y  +  part->v  -  cel->y[fc]) 

*  (cel->x[IRING(fc-f-l)]  -cel->x[fc])  -(part->x  +  part->u 

-  cel->x[fc])  *(cel->y[IRING(fc+l)]  -cel->y[fc])) 

/  ((cel->y[IRING(fc+2)]  -cel->y[fc])*(cel->x[IRING(fcTl)] 

-  cel->x[fc])  -  (cel->x[IRING(fc-h2)]  -cel->x[fc]) 

*  (cel— >y[IRING(fc-hl)]  — cel— >y[fc]))  >—  0.0)  continue;  70 

if  (tjnt  <  0)  tjnt  —  0; 

If  this  is  the  minimum  time  so  far,  record  it  and  current  face 
if  (tjnt  <  t_min){ 
t_min  =  tjnt; 
face  =  fc; 

} 

} 

I  ****  If  the  time  to  collide  with  all  faces  is  greater  than  time  remaining 
for  movement,  displace  particle  along  trajectory  and  mark  it  as  maintaining 
its  current  cell  ****/  so 
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if  {t_min  >—  time_intervai){ 

part->x  +=  part->u  *  timejnterval; 
part->y  +=  part->v  *  timejnterval; 
part->dest_cell  —  1; 
part->mvmt_time  =  0.; 

} 

!****  If  the  min  time  for  collision  is  smaller  than  remaining  time,  particle 
has  hit  something,  determine  what  that  was  and  act  appropriately  ****! 
else{ 

if  (cel~>nbr[face]  SOLID.BOUNDARY){  90 

f  If  a  solid  boundary  was  contacted,  displace  particle  to  boundary 
and  call  function  to  process  reflection  ****! 
part“>x  +=  part->u  *  t_min; 
part->y  part->v  *  t_min; 

part->mvmt_time  —  timejnterval  —  t_min; 
part->crosdJace  =  (char)face; 

Solid_Boundary(part,  cel,  face); 

} 

else  if  (cel->nbr[face]  ==  INFLOW.OUTFLOW){ 

#ifdefIO.BOUNDARY  loo 

f****  If  particle  left  grid,  mark  its  place  as  empty  and  decrement 
the  number  of  particles  in  this  cell  ****/ 
part— >dest_cell  —  0; 
cel— >n_particles — ; 
cel — >  n_vacanciesH- + ; 

#else 

1**"^*  If  we’re  not  doing  inflow/ outflow,  make  these  bounds  solid  ****/ 
part— >x  +—  part— >u  *  t_min; 
part->y  +—  part->v  *  t_min; 

part“>mvmtjime  =  timejnterval  -  t_min;  no 

part->crosdJace  =  (char)face; 

Solid_Boundary(part,  cel,  face); 

#endif 

} 

else{ 

/  *  If  it  didn’t  hit  a  boundary  of  some  sort,  particle  has  crossed 
into  another  cell  */ 

/****  Displace  particle  to  boundary  ****/ 
part— >x  +=  part— >u  *  t_min; 

part— >y  +=  part— >v  *  t_min;  120 

/****  Increment  the  current  cell  traveler  counter  ****/ 
cel->n_tvlrs++; 

/****  Decrement  the  current  cell  particle  counter  ****/ 
cel  -  >  n_part ides  — ; 

/****  Find  current  cell  index  ****/ 
curr_cell  =  cel  —  cell; 

j  Fmd  index  of  intersection  face  in  new  cell  ****/ 
for(fc  =  0;  fc  <  3;  fc++) 

if(cell[cel->nbr[face]].nbr[fc]  ==  curr_cell)  break; 
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I  ****  Store  index  of  crossed  face  in  new  cell  ****/ 
part->crosd_face  —  (char)fc; 

!****  Store  destination  cell 
part->dest_ceil  =  cel ->nbr  [face]; 

/^***  Store  time  remaining  for  movement  ****/ 
part  — >mvmt_time  =  timejnterval  -  t_min; 


Appendix  B 

Particle  Communication  Function 


This  appendix  contains  the  function  which  transfers  particles  between  cells.  It  is 
intended  to  serve  only  as  an  annotated  listing;  the  full  description  of  intercelb  com¬ 
munication  is  given  in  Section  2.6. 

Arguments 

The  communication  function  takes  two  arguments.  The  first  argument,  ‘cel’  is  a 
pointer  to  the  cell  that  the  particle  is  departing.  The  second  argument,  ‘pJndx’,  is 
the  index  of  this  particle  on  the  list  in  its  current  cell.  It  may  be  noted  that  the 
function  accesses  other  cells  in  the  course  of  its  task.  This  is  possible  because  the  cell 
data  was  made  globally- accessible  to  shorten  the  argument  lists  of  oft-called  functions. 

Code  Listing 

void  Communicate(struct  cell_unit  *cel, unsigned  short  pJndx) 

unsigned  short  d_cell,  dcpjndx; 
struct  atom  *dc_part; 

j  ****  Decrement  the  traveler  counter  for  old  cell 
cel— >n_tvlrs — ; 

I  ****  Note  this  particle's  intended  destination  ****/ 
d_cell  =  cel ->  par  tide  [pJndx].  dest_cell; 

/  ****  Mark  destination  cell  as  having  received  a  traveler  ****/ 
cell[d_cell].tvlr_recd  =  1; 


/  ****  Mark  this  particle  as  stationary  so  it  is  not  moved  by  a  subsequent  recursive  lo 

call,  leading  to  an  infinite  loop  ****/ 
cel— >  par  tide  [pjndx].dest_cell  =  1; 

!  ****  Don't  bother  searching  through  particle  list  if  this  cell  has  no 
vacancies  and  no  travelers  ****/ 
if  (cell[d_cell].n_vacancies  !=  0  ||  cell[d_cell].n_tvlrs  !=  0){ 

I****  Cycle  through  destination  cell's  particles,  looking  for  an  open  space  ***/ 


67 


for  (dc_part  —  cell[d_ceil]. particle;  dc_part  <  cell [d_cell]. particle 
+  cell[d_cell].plistjength;  dc_part++){ 

!****  If  destination  cell  for  a  particle  is  one,  its  space  is  unavailable 
(it’s  not  leaving  current  cell)  so  keep  looking  ****( 
if  (dc_part->dest_cell  — =  1)  continue; 

//  dest.  cell  is  zero,  space  is  open.  Place  particle  and  return  ****) 
else  if(dc_part— >dest_cell  ==  0){ 

Place  particle  in  new  cell  ****( 

*dc_part  =  cel->particle[pjndx]; 

Increment  new  cell’s  particle  counter  ****) 
cell[d_cell].n_particlesH-+; 

Decrement  new  cell’s  vacancy  counter  ****) 
cell  [d_cell]  .n_vacancies — ; 

Mark  its  old  position  as  vacant  ****/ 
cel— >  particle  [p_mdx].dest_cell  =  0; 
cel  -  >  n_vacancies-j-  -f ; 

return; 

} 

!****  If  dest.  cell  is  a  true  cell,  recursively  call  this  function  to  move 
the  particle  to  its  new  cell  then  put  current  particle  in  its  place  ****/ 
else{ 

dcpjndx  =  dc_part  -  cell[d_cell]. particle; 

Communicate(dcp_indx,  &cell[d_cell]); 

cell[d_cell].particle[dcpjndx]  ==  cel ->  particle  [pjndx]; 

cell[d_cell].n_particles+H-; 

cell[d_cell]  .n_  vacancies — ; 

cel->particle[p_indx].dest_cell  -  0; 

cel  -  >  n_vacancies+ + ; 

return; 

} 

} 

} 

/  *  If  no  spaces  are  found  in  list,  place  particle  at  the  end 

!****  First  ensure  there  is  room,  making  some  if  necessary  ****/ 
if  (cell[d_cell].plist_length  =—  cell[d_cell].pll_max){ 
cell[d_cell].pll_max  +=  P_LIST_INCREMENT; 
cell  [d_cell].  par  tide  =  (struct  atom  *)realloc  (cell  [d_cell]. particle, 
(size_t)(cell[d_cell].pll_max  *  sizeof(struct  atom))); 

} 

if(cell[d_cell]. particle  NULL)  printf(" FAILED  REALLOC ! \n'‘);fflush(NULL); 

/  ****  ^ow  place  particle,  increment  the  particle  counter  of  destination  cell 
and  decrement  traveler  counter  of  old  cell  ****! 
celi[d_cell]  .particle[cell[d_cell]  .plist  Jength+-h]  =  cel-  >  particle  [pjndx] ; 
cell  [d_cell] .  n_par  tides  H—f ; 
cel“>particle[pjndx].dest_cell  =  0; 
cel  -  >n_vacancies++; 
return; 


Appendix  C 

Inflow/Outflow  Function 


This  appendix  contains  the  function  which  enforces  boundary  conditions  at  inflow  / 
outflow  cells.  It  is  intended  to  serve  only  as  an  annotated  listing;  the  full  description 
of  the  inflow  outflow  boundary  treatment  is  given  in  Section  2.7.2. 

Arguments 

The  10  boundary  function  takes  only  one  argument.  This  is  a  pointer  to  an  array  of 
structures  containing  10  cells,  named  ‘io-cell’.  Each  of  these  structures  contains  the 
computed  mean  velocities  as  well  as  a  pointer  to  the  cell  in  the  main  array  to  which 
the  10  face  belongs. 

Code  Listing 

void  Process  JO  (struct  io  *io_cell) 

{ 

float  u,  u_mean,  u_norm,  u_ext,  u_distmax,  dist_max,  dist,  xtemproot, 
xjemp,  yjemp,  zjemp,  temp,  v_mean,  w_mean; 
int  inp_diff,  sgn,  flag,  pl_pos,  ip; 
unsigned  long  which; 
struct  io  *io_cl; 
struct  cell_unit  *cel; 
struct  atom  *part; 

j  ****  Enforce  boundary  conditions  at  10  cells  ****/  lo 

for  (io_cl  =  io_cell;  io_cl  <  io_cell  4-  njocells;  io_cl++)  { 

I****  Figure  out  to  which  cell  this  10  face  corresponds  **"^*1 
cel  =  &cell[io_cl->cell]; 

u_mean  —  io_cl->u_mean; 

I  ****  Determine  sign  of  inward— facing  normal 
if  (cel->x[l]  <  0.5*x_max){ 

/  ****  If  this  is  an  inflow  face,  enforce  temperature  and  transverse  speed  ****/ 
sgn  =  1; 

xjemp  1=  yjemp  =  zjemp  =  1.0;  20 

v_mean  =  w_mean  =  0; 
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temp  ~  1: 


} 

else{ 

I  ****  If  this  case  is  exhomting  to  vacuum,  skip  outflow  faces 
#ifdefVACOUT 

continue; 

#endif 

/****  If  this  is  an  outflow  face,  get  temp  and  t.v.  speed  from  flow  sample  taken 

elsewhere  (after  the  collision  routine)  ****/  30 

sgn  -1; 

v_mean  =  io_cl— >v_mean; 

wjnean  =  io_cl“>w_mean; 

xjemp  =  2*(io_cl— >u2  —  u_mean*u_mean); 

yjemp  =  2*(io_cl“>v2  -  v_mean*v_mean); 

z_temp  =  2*(io_cl->w2  -  w_mean*w_mean); 

temp  =  (x_temp  +  yjemp  -(-  zjemp)  /  3.0; 

if  (temp  ==  0.0)  temp  =  1; 

} 

I****  Calculate  difference  between  target  and  actual  number  of  particles,  40 

adding  any  round-  off  from  previous  steps  and  storing  that  from  this 
step  to  lessen  its  effect  in  overall  number  of  particles  introduced  **/ 
io_cl— >fnp_diff  H-=  cel— >n_particies  - 
(io_cl“>pressure  *  cel->volume  *  n_inf  /  temp); 
inp_difF  =  io_cl->fnp_diff; 
io^cl— >fnp_difF  —  =  inp_diff; 

/  *  Add  particles  if  there  aren’t  enough  in  cell  */ 

if  (inp_diff  <  0)  { 

!  ****  Normal  velocity  is  positive  for  inflow  ****/ 

u_norm  =  sgn  *  u_mean;  50 

1*^**  Calculate  the  maximum  of  the  incoming  velocity  distribution’^***/ 
u_distmax  =  sgn*0.5*(— u_norm  4-sqrt(u_norm*u_norm  +2)); 
dist_max  =  sgn*  (u_distmax+sgn*u_norm)  *exp  ( — (u_distmax)  *  (u_distmax)) ; 

xtemproot  sqrt(x_temp); 

/****  Set  the  extremum  of  particle  speed.  If  u  and  ujiorm  are  equal  and 
opposite,  the  particle  will  not  cross  the  boundary  (inflow  only)**/ 
if  (u_norm  <  3.0) 

u_ext  —  —  u_norm; 

else 

u_ext  =  — 3.0*xtemproot;  /* restrict  the  selection  range  to  lower  rejections*/  60 

/****  First,  ensure  there  is  room  on  cell’s  list  for  incoming  particles  ***/ 
if  (cel->plist_length  H — inp_diff  —  cel->n_vacancies  >=  cel->pll_max){ 
cel->pll_max  +=  2  *  (-inp_diff  -cel->n_vacancies); 
cel— >particle  =  (struct  atom  *)realloc(cel— >particle, 

(size_t)(cel->pll_max  *  sizeof(struct  atom))); 

} 

/****  Initialize  particle  list  position  pointer  ^^**/ 
pLpos  =  0; 
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for  (ip  —  0;  ip  <  -inp_diff;  ip++){ 

Search  remaining  positions  in  particle  list  for  empty  spots  ****/ 
for(  ;  pLpos  <  cel— >plist_length;  pl_pos++) 
if  (cel->particle[pl_pos].dest_cell  ===  0){ 
cel  -  >n_vacancies — ; 

break; 

} 

!  ****  Incerement  the  cell’s  particle  counter  ****/ 
cel  -  >  n_part  ides + H- ; 

I  ****  If  above  search  put  us  off  particle  list,  extend  it 
if  (pLpos  >=  cel->plist_length)  cel->plist_length-|-H-; 

I  Assign  a  pointer  to  this  particle  for  later  convenience  ****!  so 

part  =  &:cel->particie[pl_pos]; 

Increment  particle  list  position  for  next  cycle  ****! 
pLpos++; 

Set  particle’s  destination  cell  and  crossed  face  registers 
part->dest_cell  =  1; 
part->crosd_face  =  ( char  )io_cl->  face; 

Place  the  particle  randomly  on  the  cell’s  10  edge  ^***1 
part->x  =  cel- >x[io_cl->  face]; 
part->y  =  cel- >y[io_cl->  face]  +  ran2(seed)* 

(cel— >y[IRING(io_cl— >  face  +  1)]  —  cel— >y[io_cl— >  face]);  90 

!****  Assign  tangential  velocities  according  to  equilibrium 
distribution  ****/ 

part->v  ~  v_mean  +  sqrt(y_temp*-log(ran2(seed)))  *  sin(2.0  *  M_PI 
*ran2(seed)); 

part->w  =  w_mean  +  sqrt(z_temp*— log(ran2(seed)))  *  sin(2.0  *  M_PI 
*ran2(seed)); 

j****  Set  u  using  a  selection/ rejection  technique  on  a  fiuxal 
distribution  ****/ 

do 

{  100 
/  ****  Randomly  select  a  value  of  u  ****/ 
if  (sgn>0) 

u  =  (u_ext  +  (3.0  *xtemproot  -  u_ext)  *ran2(seed)); 

else 

u  =  (-(u_mean  +  3.0  *xtemproot)  +  3.0  *xtemproot  *ran2(seed)); 

/  ****  Compute  the  value  of  the  distribution  for  this  u  ****/ 
dist  =  sgn*(u/xtemproot  -h  u_mean)*exp(— (u*u)/x_temp)/dist_max; 

jwhile  (dist  <  ran2(seed)); 

/  *^**  Add  mean  velocity  to  assigned  thermal  speed 

part— >u  =  u_mean  +  u/xtemproot;  no 

/  ****  Move  new  particle  a  random  fraction  of  one  time  step  ****/ 
part->mvmtjime  =  ran2(seed)  *  T_STEP; 

Move_Particle(cel,  part,  1,  SELECTIVE_MOVE); 
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} 

} 

Process  comm, unication  links  until  all  are  empty 

do{ 

flag  0; 

for(cel  =r  &cell[2];  cel  <  &cell[2]+N_CELLS;  cel-f+){  12c 

!****  Skip  if  this  cell  has  no  travelers  ****/ 
if  (cel->n_tvlrs  ==  0)  continue; 

Set  a  flag  to  signify  that  a  traveler  was  found  ****/ 
flag  =  1; 

Process  all  travelers  in  this  cell  ****! 

/*NOTE:  Do  not  use  pointers  into  the  particle  array  in  Communicate 
or  its  calls  because  it  contains  a  realloc* / 
for(ip  =  0;  ip  <  cel->plist_length;  ip++){ 

j  ****  Skip  if  this  is  an  empty  position  or  stationary  particle 
if  (cel”>particle[ip].dest_cell  <  2)  continue;  iso 

Communicate(ip,  cel); 

I  ****  Go  on  to  the  next  cell  if  we  just  placed  the  last  traveler  ****/ 
if  (cel->n_tvlrs  ==  0)  break; 

} 

} 

j****  If  communication  was  performed,  move  newly— placed  travelers  ****j 
if  (flag){ 

for(cel  =  &cell[2];  cel  <  &cell[2]+N_CELLS;  cel++){ 

j****  Skip  this  cell  if  it  received  no  travelers  ****/  140 

if  (!cel— >tvlr_recd)  continue; 

j****  Reset  cell  ’traveler  received’  flag  ****f 
cel— >tvlr_recd  =  0; 

/  ****  Move  any  particles  which  have  remaining  time 
Move_Particle(cel,  cel->particle,  cel->plist_length,SELECTIVE_MOVE); 

} 

}while(flag); 

return; 

}  150 
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